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Then be it ours with steady mind to clasp 
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Chapter 1 

General Introduction 



Nature has a hierarchical structure, with time, length, and energy scales rang- 
ing from the microscopic to the macroscopic. Surprisingly, it is often possible 
to discuss these levels independently, e.g., in the case of molecular fluids. The 
macroscopic level, i.e. the world directly perceived at our human senses, is 
described hydrodynamically. That is, by continuous (or piecewise continuous) 
functions of the spatial coordinates r and time t (hydrodynamic fields). As 
a result, the great disciplines of macroscopic physics, such as fluid mechanics, 
elasticity, acoustics, electromagnetism, are field theories. These fields are deter- 
mined by integro-differential equations involving unknown functions of r and 
t. The microscopic level is described as a collection of a very large number of 
constituents, such as atoms, molecules, or other more complex entities inter- 
acting with each other through well-defined electro-mechanical forcefl Their 
evolution with time is provided by the laws of quantum mechanics, but in many 
cases classical mechanics results a good approximation. 

If there is a clear scale separation between microscopic and macroscopic, 
then an intermediate scale can be introduced. Such a level of description is 
referred to as the mesoscopic level. It is a coarse-grained representation to cap- 
ture the physics of microscopic models on the large length and time scales, which 
is described by stochastic partial differential equations [HHI. By concreteness, 
the mesoscopic description of the system in question, whose dynamics is on the 
average governed by the laws of macroscopic time evolution, is a result of build- 
ing up of microscopic fluctuations. This approach aims at understanding the 
slowly- varying long wavelength, low frequency properties of many-body systems 
and form the basis of the field theoretical analyses of critical phenomena. In the 
literature it is common include in this intermediate scale the kinetic description 
of gases, e.g., Boltzmann equation. Thus, this level is also referred to as the 



^The difference between microscopic and macroscopic levels is essentially relative. The key 
concept here is the large number of constituents, not their small size. For instance, a galaxy 
(macroscopics) is a object composed by a large number of stars (microscopies), a colony of 
organisms is a large set of individuals, a gas comprised by a large number of molecules, etc. 
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kinetic levejl. This is in fact a reduced microscopic level, described by one- 
particle coarse-grained distribution functions and the asymptotic, irreversible 
kinetic equations. 

The three different levels of descriptiord are, in principle, equivalent ap- 
proaches of the same reality. Nevertheless, the form of the laws derived at these 
three levels are so different from each other that there soon appeared a stringent 
need for an explanation of hydrodynamics in terms of the microscopic evolu- 
tion of the underlying collection of constituents. The bridges between them are 
provided by Statistical Mechanics. Statistical mechanics aims to derive macro- 
scopic behavior of matter originated in the cooperative behavior of interacting 
(microscopic) individual entities. Some of the phenomena are simple synergetic 
effects of the actions of individuals, e.g., the pressure exerted by a molecular 
gas on the walls of its container and thousands of fireflies flashing in synchrony; 
while others are paradigms of emergent behavior, having no direct counterpart 
in the properties or dynamics of individual constituents, e.g., the transition from 
laminar to turbulent flows in fiuids and the one million of atoms which give rise 
to the program of life: the deoxyribonucleic acid or DNA for short. In the latter 
cases, the behavior of the constituents become singular, but very different from 
what they would exhibit in the absence of the others. 

The most successful achievement of statistical mechanics is Ensemble Theory 
[3], which yields the connection between the macroscopic properties of equilib- 
rium systems from the laws governing the microscopic interactions of the in- 
dividual particles. It is said that a certain system stay at equilibrium when 
it is isolated, shows no hysteresis, and reaches a steady state (all its macro- 
scopic properties are fixed) [1]. In that case, the macroscopic properties are 
expressed in terms of thermodynamic variables. Indeed, the thermodynamics is 
in this sense a kind of hydrodynamic description, which consists only of laws and 
relations between thermodynamic quantities. From the mathematical point of 
view, ensemble theory can be completely axiomatized. This provides us, at least 
in principle, analytic expressions for these thermodynamic quantities allowing 
us to derive all the relevant macroscopic information of the thermal system in 
questior0, e.g., the free energy, equation of state, etc. 

In nature, by contrast, equilibrium systems are an exception (even an ideal- 
ization) rather than the rule: non- equilibrium phenomena are overwhelmingly 
more abundant. Galaxies, human beings, chemical reactions, geophysical fiows, 
carriers in semiconductor devices, financial stock markets, ratchet-effect trans- 
port, traffic flow on a highway, to cite just a few, are (many-body) systems 



■^Formally, both stochastic differential equations and kinetic equations may belong to the 
same level, although they involve different (mesoscopic) length and time scales. 

''There may be larger scale structures (e.g., the Karman vortex train) produced by a large 
scale fluid dynamics, but we do not pay attention to such large scales. 

*Once the microscopic Hamiltonian (say T-L) of the system is specified the "canonical" sta- 
tionary probability distribution V is known in terms of the Boltzmann factor V = e^^/'^B^/Z, 
where Z is the partition function, kg is the Boltzmann's constant, and T the temperature. 
Thus, averages over such a distribution of time-independent observables can be computed. 
The remaining difficulties are merely technical. 
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Figure 1.1: Complexity arises at all levels. Some amazing examples are the 
following. Upper left: river networks in Yemen (source NASA). Upper right: 
aeolian patterns in Colorado, USA (photo by Bob Bauer). Lower left: tornado 
as a rotating geophysical flow. Lower right: oscillons in a colloidal suspension[S]. 



under nonequilibrium condition^0. Out-of-equilibrium systems are character- 
ized by not being closed, i.e., by having an exchange of energy, particles, and/or 
information with their environment. In general, the state of a nonequilibrium 
system is not determined solely by external constraints, but also depends upon 
its history (hysteretic systems). Consequently, this gives rise to the enormous 
complexity of our world which occurs to all levels of description [5] and man- 
ifests in pattern formation, absorbing states, self-organization, chaos, ageing, 
avalanches, morphogenesis, oscillations, fractals, etc [3 [HI HI [TO]. All these 
amazing and complicated phenomena are associated to instabilitiei £|, which are 
variously described as nonequilibrium phase transitions [7j llli I12i 113] , bifurca- 
tions or synergetics, with the aim of connecting microscopies with the coherent 
structures observed. 

However, most of studies of nonequilibrium systems adopt a macroscopic 
(phenomenological) point of view. In fact, little is known in general on why are 
there such an interesting structures or how this complexity arrises from interac- 
tion at a microscopic level. Consequently, no theory exists and the development 

^Nonequilibrium phenomena are also encountered whenever systems are relaxing towards 
an equilibrium steady state. Nevertheless, in most of cases these can be understood in terms 
of equilibrium concepts. 

'^Nonequilibrium instabilities are attended by ordering phenomena analogous to those of 
equilibrium statistical mechanics; one may therefore refer to nonequilibrium phase transitions 
Although nonequilibrium phase transitions represent a more varied picture than their 
equilibrium counterparts. 
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of a solid theoretical background is nowadays in an early stage compared with 
equilibrium, where the ensemble theory holds successfully. As a matter of fact, 
the key point of difference between equilibrium and nonequilibrium statistical 
mechanics is that whereas in the former case the stationary probability distribu- 
tion is known, out of equilibrium one must find the time-dependent distribution, 
which obey a general evolution equation, i.e., a master equation. In most of sol- 
uble cases — which in fact are just a few — , this can be done only approximately. 
One can only provide some unified view when the systems are not too /otQ from 
equilibrium |14| . where the system can be treated perturbatively around the 
equilibrium state and studied using linear response theory. Nevertheless, our 
attention here lies at systems far from equilibrium where such schemes break 
down. 

Within the framework given above, the main subject of this thesis 
rests on the study — at different levels of description — of instabilities 
in systems which are driven, i.e., maintained far from equilibrium 
by an external forcing. We focus here on two main classes, namely, 
driven diffusive fluids and driven granular gases. 

The first family [driven- diffusive fluid^ corresponds to systems which ar^ 
coupled to two reservoirs of energy in such a way that there is a steady energy 
flow through the system. Clearly, this definition comprises a bunch of very var- 
ied systems. We restrict ourselves to systems in which a (non-zero) current of 
particles (whose number is a conserved quantity) is set up through the system, 
there is spatial anisotropics associated with an external driving field, and even- 
tually a nonequilibrium steady state is reached. The simplest example seems to 
be a resistor gaining energy from a battery and losing it to the atmosphere. But 
even for this reduced class of systems the stationary probability distribution is 
unknown. 

The basic motivation behind the introduction of driven- diffusive fluids cor- 
responds to a need of unravel the key ingredients of nonequilibrium behavior. 
A reasonable approach consist in investigating systems as simple as possible, 
which aim for capture the microscopic essentials yielding to the complicate, 
macroscopic nonequilibrium ordering. These microscopic models are usually 
oversimplified representations of real systems, and consider entities as particles 
interacting via simple rules. Between them, lattice models — which are based in 
discretization of space into lattice sites — have played a important role due to 
the fact that they sometimes allow for exact results, and allow one to isolate the 
specific features of a system. In addition, lattice models allow us to obtain the 



^Nevertheless, a precise definition of not too far is rather difficult to give; the borderline 
is rather unclear. 

*We should mention here that in the literature these systems are commonly referred as 
driven— diffusive systems instead of driven-difFusive fluids, which is adopted in this thesis. 
As we shall describe in some detail in Chapter [3] we reserve the former denomination for a 
mean-field mesoscopic description. 

^Following Schmittmann & Zia in Ref. |15| . 
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intuition one needs in order to develop the appropriate theoretical tools and, 
from a experimental point of view, are easier to be implemented in a computer. 
Such an approach has given rise to a lively research activity in the last decade 
pn [TCI [TBI [T71 [T51 [TOl [20] and a bunch of emerging techniques may now be 
applied to lattice systems, including nonequilibrium statistical field theory. A 
general amazing result from these studies is that lattice models often capture 
some essentials of social organisms [HI [HI [13] , formation of river networks [H] , 
epidemics [53], glasses, electrical circuits, traffic [551 HI], hydrodynamics, col- 
loids and foams, enzyme biology [5H], living organisms [Hj, and markets [50] . 
for example. In particular, lattice models have succeed in understanding (equi- 
librium) phase transitions and critical phenomena [16l [31] [32] due, in many 
aspects, to renormalization group methods [33l 134) . 

The central result in (equilibrium) critical phenomena is universality, i.e., 
the behavior of disparate systems in the vicinity of critical points are deter- 
mined solely by basic features — dimensionality, range of forces, symmetries, 
etc. — and does not seem to depend to any great extent on the particular sys- 
tem [33l [35l [36] . One can therefore hope to assign all systems to classes each 
of them being identified by a set of critical indices (exponents and some am- 
plitude ratios). How much of these concepts apply to nonequilibrium phase 
transitions are only just emerging, and one expects lattice models to be equally 
important here. These issues are considered in Chapter [3J in which a particular 
driven-diffusive lattice model, prototype for nonequilibrium phase transitions, 
is investigated. However, a well-known disadvantage of lattice models is that 
when they are compared directly with experiment result too crude. This has 
to be understood in the sense that they often do not account for important 
features of the corresponding nonequilibrium phase diagram, such as structural, 
morphological, and even critical properties. Furthermore, theoreticians often 
tend to consider them as prototypical models for certain behavior, a fact which 
is in many cases not justified. This will be discussed in Chapter [4] where we 
introduce a novel, realistic model for computer simulation of anisotropic fluids. 

The second class of systems we consider in this thesis concerns driven granu- 
lar gases. A granular material is a large conglomerate of discrete macroscopic^ 
(classical) particles. Examples, which occur at very different length scales, may 
range from powders to intergalactic dust clouds; they include sand, concrete, 
rice, volcanic flows, Saturn ring's, and many others (for a review, see [37]). They 
have important applications in industrial sectors as pharmaceutical, construc- 
tion and civil engineering, chemical, food and agricultural, etc. On a larger 
scale, they are also relevant in understanding many geological and astrophysical 
processes. In addition, they are closely related to a broader class of systems, 
such as foams, colloids and glasses. The collective behavior of an assembly of 
granular particles exhibits an impressive variety of phenomenelll] which include 
a plethora of pattern forming instabilities [38l [Ml EHl SI] , clustering [42l l43l [44] . 

^''In the sense that they are directly perceived at our human eyes. 

^^The size range of phenomena goes from lO^^ra in powders, to ~ 10'^ m in deserts, up 
to lO^^m in planetary rings. 
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flows and jamming in hoppers [45], avalanches and slides [46], mixing and seg- 
regation [47], convection and heaping [48l [49], eruptions [50] [51], to cite just 
a few. Clearly, their relevancy is beyond any doubt, and understanding their 
properties is not only an urgent industrial need, but also an important challenge 
for physicists. 

Since grains have a macroscopic size, friction and restitutional losses from 
collisions give rise to dissipative interactions — kinetic energy is continually 
transferred into heat — . Although the dominant interactions depend on particle 
size. The relevant energy scale for a typical, say, rice grain of length I « 5mm 
and mass m is its potential energy mgl, where g is the gravitational accel- 
eration. For much smaller grain sizes other kind of interactions may become 
importantly. Hence, in granular materials the thermal energy ksT is irrele- 
vant, i.e. kfiT ^ mgl. This implies that dynamical effects outweigh entropy 
considerations [37 and, therefore, neither equilibrium statistical-mechanics nor 
thermodynamics arguments are useful. Furthermore, due to the fact that gran- 
ulates are comprised by a finitenumber of particleqlfjusually there is a strong 
effect of fluctuations and, therefore, statistics may be dominated by rare events. 
Other important effects involve the interstitial fluid, e.g. air or water, although 
in many situations one can ignore it with confldence and consider the grain- 
grain interactions alone. In this thesis (specifically in Part [IT]) we will restrict 
ourselves to dry granular systems, where the interaction between particles take 
place mainly via contact forces. 

However, in spite of the simple (classical) features which characterize them, 
granular materials behave in a rather unconventional way: their phases — solid, 
liquid, and gas — have complex-collective properties that distinguish them from 
molecular fluids and solids. As a result, their statistics and dynamics are poorly 
understood. In fact, from a theoretical point of view, these systems are still 
only understood in terms of predictions of a general nature and many open 
questions remain. Linear elasticity does not apply for granular solids, which are 
highly histeretic, show static stress indeterminacy, and possess static yield shear 
stress. The interlocking of grains leads to force chains, jamming, dilation on 
shearing, etc. Liquids are viscous, show avalanches (possess a critical slope), size 
segregation in mixtures, shear bands (narrow regions separating blocks moving 
with different velocities), pattern formation. Gases, which can only persist 
with continuous energy supply, are highly compressible, show clustering, long 
range correlations, inelastic collapse in computational models, non-maxwellian 
distribution functions, pattern formation, and lack of scale separation between 
macroscopic scales and microscopic ones. All of this phenomenology makes 
difficult to cast them in the framework of the classical three different phases of 
matter. Some of the interesting open issues and controversies that we address in 
Chapters [6| and [7| are: What are the statistical properties of granular systems. 



^■^For instance, it may occur charging and surface coating for d ~ 0.1, or magnetization and 
surface adhesion for d ~ 0.01. 

^•^In contrast with molecular systems N << N^, where is the Avogadro number. Nev- 
ertheless, in certain cases a few thousands or even hundreds of particles are enough to allow 
for statistical approaches. 



1.1 Overview 
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how do granular phase changes occur, and what are the optimum continuum 
models and when do they apply. In particular, we focus on the continuum 
description of clustering, symmetry breaking, and phase separation instabilities 
in granular gases. As we shall see, these instabilities provide sensitive tests 
to models of granular flows and contribute to the understanding of pattern 
formation far from equilibrium. 

1.1 Overview 

In this thesis, we investigate the two aforementioned types of systems, which 
have become of paramount importance in recent years. The study is divided 
into four parts. The Part HI which comprises Chapters [51 [3] and SI is devoted to 
driven- diffusive fluids. Part HIl deals with driven granular gases, and consists of 
Chapters [3 H and H Part HII] includes the Appendix El [B O and Finally, 
we close (Part ITV)) with a detailed summary and the main conclusions. 

The first part begin in Chapter [2]bv introducing the driven lattice gas — 
prototype for nonequilibrium anisotropic phase transitions — . This model may 
be viewed as an oversimplification of certain traffic and flow problems. The par- 
ticles interact via a local anisotropic rule, which induces a preferential hopping 
along one direction, so that a net current sets in if allowed by boundary condi- 
tions. We discuss some of its already known properties and controversies. We 
shall describe in greater detail the computational and field-theoretical methods 
employed in Chapters [3] and [H 

One of the most important questions we can ask about any model is whether 
the behavior that it displays is universal. In Chapter [3]we address this question 
and study the phase segregation process in the prototypical lattice (microscopic) 
model introduced in Chapter [51 The emphasis in our study is on the infiuence 
of dynamic details on the resulting (non-equilibrium) steady state, although we 
also pay some attention to kinetic aspects. In particular, we shall discuss on 
the similarities and differences between the driven lattice gas and related lattice 
and off-lattice models, from a simulational and theoretical point of view. The 
comparison between them allows us to discuss some exceptional, hardly realistic 
features of the original driven lattice gas. In addition, we test the validity of 
the two competing mesoscopic (Langevin-type) approaches on describing the 
critical behavior of these models. 

We introduce in Chapter [4l a novel, realistic off-lattice driven fluid for 
anisotropic behavior. We describe short-time kinetic and steady-state prop- 
erties of the non-equilibrium phases, namely, solid, liquid and gas anisotropic 
phases. This is a computationally-convenient model which exhibits a net cur- 
rent and striped structures at low temperature, thus resembling many situations 
in nature. We here focus on both critical behavior and some details of the nu- 
cleation process. Concerning the critical behavior, we discuss on the role of 
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symmetries on computer and field-theoretical modeling of non-equilibrium flu- 
ids. In addition, some important comparisons with experiments are drawn. 

Chapter [5] serves as an introduction to PartUll which is devoted to driven 
granular gases. We briefly outline the current situation with regard to research 
on granular gases. In particular, we discuss on the applicability of a hydrody- 
namic theory to granular flows. We also present the basic ingredients of the 
studied models and the computer simulation methods. 

In Chapter [6] we address phase separation instability of a monodisperse 
gas of inelastically colliding hard disks confined in a two-dimensional annulus, 
the inner circle of which represents a "thermal wall" . Our main objectives are 
to characterize the instability and compute the phase diagram by using granu- 
lar hydrodynamics — we employed for hydrodynamics Enskog-type constitutive 
relations, which are derived from first principles — and event driven molecular 
dynamics simulations. We also discuss on clustering and symmetry-breaking 
instabilities. By focusing on the annular geometry, we hope to motivate experi- 
mental studies of the granular phase separation which may be advantageous in 
this geometry. 

This Chapter [7] addresses granular hydrodynamics and fiuctuations in a 
simple two-dimensional granular system under conditions when existing hydro- 
dynamic descriptions break down because of large density, not large inelasticity. 
We study a model system of inelastically colliding hard disks inside a circular 
box, driven by a thermal wall in which a dense granular cluster behave like a 
macro-particle. Some features of the macro-particle are well described by the 
stationary solution of granular hydrostatic equations by using phenomenologi- 
cal constitutive relations. Hydrostatic predictions arc tested by comparing with 
molecular dynamics simulations. We are able to develop an effective Langevin 
description for the macro-particle, confined by a harmonic potential and driven 
by a delta-correlated noise. We are also concerned about whether the crystal 
structure depends on the geometry. 

The Part IIIII is comprised by four appendices: Appendix [A] is devoted to 
the most important properties of the Lattice Gai^. which is the equilibrium 
counterpart of the driven lattice gas studied in Chapter |3l in Appendix [B] 
we introduce a model which aim at characterizing the stability of the interface 
in driven lattice models; Appendix [C] details some field theoretical deriva- 
tions performed in Chapter [3] concerning Langevin-type equations; and in Ap- 
pendix [D] we propose a set of constitutive relations for the hydrostatic descrip- 
tion of three-dimensional granular gases. 

Finally, after the appendices, we present in Part IIVI our main conclusions. 



^■'The lattice gas is the conserved version of the Ising model | 32l I52| . As a result, they are 
isomorphic. 
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summing up our original contributions and pointing out the future work that 
should be addressed. Also we provide a list of publications associated with the 
work performed in this thesis. 
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Part I 

Driven Diffusive Fluids 



Chapter 2 

Introduction 



As outlined in the general introduction (Chapter [1} , the first part of this thesis 
is devoted to the study of nonequilibrium phase transitions in Driven Diffusive 
Fluids. This is a broad class of nonequilibrium systems which are characterized 
to be coupled to two reservoirs of energy. The one injects energy continously 
whereas the other substracts it in such a way that a nonequilibrium steady state 
is reached asymptotically. Studies of these systems have generally focused on 
lattice systems, i.e., simplified models based on a discretization of space and in 
considering interacting particles that move according to simple local rules. This 
is because a general formalism analogous to equilibrium statistical mechanics is 
lacking even for nonequilibrium steady states. The way to overcome this draw- 
back is to develop and study simple models — while retaining the essence of the 
difficulties of nonequilibrium states — capable to capture the essential behavior. 

Within this context, one of most relevant models is the one devised by 
Katz, Lebowitz and spoh43 [sg. Motivated by both the theoretical interest 
in nonequilibrium steady states and the physics of fast ionic conductors [TT| . 
they introduced a minor modification to a well-known system in equilibrium 
statistical mechanics: the Ising model (321 [52]. It consists of particles — instead 
of spins — diffusing on a lattice, subject to an excluded volume constraint and 
an attractive nearest neighbor interaction. The particles, whose total number 
remains conserved, are driven by an external (constant) field E, settled on a 
thermal bath at temperature T. We shall refer henceforth to this model as the 
Driven Lattice Gas (DLG) (HKl^. The dynamics in the DLG is according to 
Metropolis [54] transition probabilities per unit time. Consequently, for peri- 
odic boundary conditions (and for any non-zero driving field), the result is a 
net current of particles that generates heat which is absorbed by the thermal 
bath. As E is increased, the system reaches saturation, i.e., particles cannot 
jump against the field. Eventually, anisotropic phase segregation occurs in the 
DLG for low enough bath temperatures: a striped liquid (rich particle) phase, 



Sometimes is referred to as the KLS model, after the initials of its inventors. 
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Figure 2.1: Typical configurations in the stationary regime of tire DLG for a 
128 X 128 half-filled (p = 1/2) lattice and (from left to right) T = 0.70, T = 0.95, 
T — 1.30, and T = 1.50. The temperatures are rescaled according with the 
equilibrium lattice gas critical temperature (see Appendix [A| : Tons = 2.269 in 
energy units. The value of the field is E = 1000, which is effectively infinite at 
all temperatures. 



which coexists with its vapor (poor particle) phase, percolates along the field 
direction. However, for high enough temperatures a disordered phase appears, 
characterized by a gas-like state. As a matter of fact, for two dimensional half- 
filled lattices the DLG undergoes a second-order phase transition; otherwise the 
transition is of first-order. Its phase behavior is best illustrated by 'snapshots' 
of typical configurations generated by Monte Carlo (MC) simulations [SS], as 
in Fig. O 

In spite of a deceptively simple definition, the steady states observed exhibit 
a rich and counterintuitive behavior, with many interesting features such as the 
violation of the fluctuation-dissipation theorem, a singular structure factor, and 
slow decay of spatial correlations. But perhaps its most striking, counterin- 
tuitive feature is that a strong field raises the critical temperature above the 
equilibrium lattice gasH critical temperature. We shall discuss these features in 
greater detail in Chapters [3] and H) The DLG has been also useful to model, for 
instance, ionic currents [11, 56^, traffic and pedestrians flows [57l[58], water- in- 
oil microemulsions [SS], electrophoresis experiments [50], etc (see Ref. [TS] and 
references therein.) Moreover, the highly anisotropic conflgurations exhibited 
by the DLG resembles actual driven systems. Many natural systems exhibit 
spatial striped patterns on macroscopic scales [H . These are often caused by 
transport of matter or charge induced by a drive which leads to heterogeneous 
ordering. Such phenomenology occurs during segregation in driven sheared sys- 
tems [611 inn [63] , in flowing fluids [64] , and during phase separation in colloidal 
[65] . granular [66l [67], and liquid-liquid [68] mixtures. Further examples are 
ripples shaped by the wind in sand deserts [69l [70] , the lanes and trails formed 
by living organisms (animals and/or pedestrians) and vehicle traffic [271 l29] . 
and the anisotropies observed in high temperature superconductors [TT] [72] and 
in two-dimensional electron gases [731 IZl] • 

As a result, this model system became the theoreticians' prototype for 
anisotropic behavior, and is nowadays the most thoroughly studied system show- 
ing an anisotropic nonequilibrium phase transition. 



^See Appendix [a] for a quick review of the lattice gas properties. 
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In Chapters [3] and [4] we report on novel findings on this and other related 
models by describing different realizations observed in (Metropolis) MC simu- 
lations. In particular, we study to what extent the DLG features are universal, 
or in other words, how robust is its behavior. The clue is that the DLG is, in a 
sense, pathological. We also devise a novel driven fluid for anisotropic phenom- 
ena, which we believe is not only of theoretical importance, but also relevant 
for experimentalists. In addition to MC techniques, field-theoretical Langevin- 
type approaches on the DLG are also developed. But before deep further on 
these matters, we shall describe in greater detail the computational and field- 
theoretical methods employed in these chapters. 

In statistical mechanics, Monte Carl(^ simulations |55| play a major role 
on the study of phase transitions and critical phenomena. Any MC algorithm 
generates stochastic trajectories in the system's phase space, in such a way 
that the properties of the system are derived from averages over the different 
trajectories. To be specific, for systems both in and far from equilibrium, MC 
simulations involve long sequences of configurations, which evolve from one c 
to another c' according to a defined transition probability per unit time (rate), 
namely, w(c — > c'). Once initial transients have decayed stationary observables 
can be computed as configurational averages. For systems under equilibrium 
conditions, such an average over configurations is referred to as the ensemble 
average. In such a case, and if the ergodic hypothesis [3. is assumed, ensemble 
averages are equivalent to time averagefQ. 

The acceptance rules for transitions between configurations are chosen such 
that these configurations occur with a frequency prescribed by the desired prob- 
ability distribution. For equilibrium systems, this must be the "canonical" sta- 
tionary Gibb's distribution Vo = e^'^^^^'^ /Z, where H is the Hamiltonian, Z is 
the partition function, ke is the Boltzmann's constant, and T the temperature. 
Since our interest here is in nonequilibrium behavior, we will need to specify 
how a given configuration c evolves into a new one, c', that is, w(c — > c'). 
As a consequence, we must deal with a time-dependent probability distribution 
function V{c,t), which obeys a master equation [HIS] 

9,p(c,t) =^{-p(c',t)w(c' ^c]-7'(c,t)w(c^c')} . (2.1) 

c' 

Its stationary solution, V(c), controls all the time-independent properties. To 
ensure that the desired equilibrium distribution Vo is reproduced, one chooses 



^More precisely, one should refer it as the Monte Carlo importance sampling algorithm 
|54l I55| . The nickname "Monte Carlo" comes from the fact that this technique entails a large 
sequence of random numbers. 

■'in fact, if MC simulations mimic configurational averages, Molecular Dynamics is a scheme 
for studying the natural time evolution (time averages) of a certain system. For thermalized 
(or equilibrated) systems differences between both techniques should disappear in the thermo- 
dynamic limit. Chapters [6] and [7] deal with Molecular Dynamic simulations, which are briefly 
described in Chapter |5] 
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rates which satisfy the detailed balance condition, namely 

Mc' ^ c] ^ Vo[c) 

w(c c'] Vo[c'] ' ^ ' ' 

This ensures that the bracket in Eq. (I2.ip vanishes. The important point here is 
that the ratio Vo[c)/'Pq[c') = e-^-H/^BJ ^ ^here AH = H[c')-n[c). One may 
thus choose rates of the form w(c — > c') — f^AH/lcBT), with an appropriate f 
satisfying 

f (-x)/f (x) = . (2.3) 

Some choices for f are the Kawasaki rate f(x) = 2(1+6") [75], the van 
Beijeren & Schulman rate f(x) — e^"^^ [76], and the Metropohs rate f(x) — 
min{1, e^"} [SI]. All of our simulations in Part U concern the Metropolis rate. 

The most simple way of driving a system (as in the DLG case) into a nonequi- 
librium steady state is to impose rates that violates detailed balance. A straight- 
forward extension of AH = T-L[c') —T-L[c) in Eq (|2.2I) is to include the work done 
by the field, i.e., to define the total energy difference of the for AH + Ex, with 
X the particle displacement in the field direction. This rate with the property 
Eq. (|2.3p satisfies the detailed balance condition locally but not globally so that 
microscopic reversibility of the process as in equilibrium is not guaranteed. This 
scheme is adopted in the DLG which (together with toroidal boundary condi- 
tions) yields the system towards a nonequilibrium steady state. 

It is worth notice that there are recent computational coarse-grained meth- 
ods, e.g., Dissipative Particle Dynamics (DPD) and Brownian Dynamics (BD) 
(see for instance [77]), for fluid dynamics which may have several advantages 
over conventional computational dynamics methods. The DPD method — first 
proposed by Hoogerbrugge and Koelman [7H] using heuristic arguments, and 
properly formalized by Espahol and Warren |79| — is a useful technique when 
studying the mesoscopic structure of complex liquids. BD is a technique that 
resembles DPD but each (mesoscopic) particle feels a random force and a drag 
force relative to a fixed backgrounq^- However, as we stressed in Chapter (U 
we are interested in how the macroscopic behavior emerges from microscopic 
(atomistic) interactions rather than in modeling the dynamics of fluids from a 
mesoscopic point of view. 

Statistical field theories are a complementary approach to the understanding 
of nonequilibrium ordering in the DLG and in related models. These approaches 
often pose great conceptual and computational challenges in themselves, but we 
will focus here on critical properties and renormalization group notions through 
the Langevin equation [T] [2] . This is a stochastic partial differential equation, 
which corresponds to a mesoscopic (coarse-grained) description of the system, 
and it is thought to describe low-frequency, small-wave number phenomena 
(such as critical properties) . It takes into account only the relevant symmetries 

^The important difference between BD and DPD is that BD does not satisfy Newton's 
third law and hence is does not conserve momentum. As a consequence, BD cannot reproduce 
hydrodynamic behavior. 



17 



for the problem under study and lets the fast degrees of freedom that we forget 
about in the coarse-grained description to act as noise. The Langevin equation 
for a field 4> takes the following form [U [2] 

^^^=F[(t)] + G[ct)]t(r,t), (2.4) 

where f, is a random variable which represents a uncorrelated Gaussian white 
noise. F and G, which are determined in general from symmetry considerations, 
are analytic functionals of 4). In fact, these functionals include every analytic 
term consistent with the symmetries of the microscopic system. In principle, 
F and G can be determined phenomenologically through experiments and/or 
symmetry arguments. Other difficult procedure which is seldom carried out 
is the coarse-graining operation that produces the Langevin equation from the 
microscopic model, which follows, say, a master equation as Eq. (|2.ip . We adopt 
the latter procedure (we extend the scheme devised in Ref. ^80^) in the Chapter[3] 
and in greater detail in the Appendix [Cl 
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Chapter 3 



Lennard— Jones and Lattice 
models of driven fluids 



The present chapter describes Monte Carlo (MC) simulations and field theoret- 
ical calculations that aim at illustrating how slight modifications of dynamics at 
the microscopic level may influence, even quantitatively, the resulting (nonequi- 
librium) steady state. With this aim, we take as a reference the driven lattice gas 
(DLG), which we described qualitatively in the previous chapter. We present 
and study related lattice and off-lattice microscopic models in which particles, 
as in the DLG, interact via a local anisotropic rule. The rule induces preferential 
hopping along one direction, so that a net current sets in if allowed by bound- 
ary conditions. In particular, we shall discuss on the similarities and differences 
between the DLG and its continuous counterpart, namely, a Lennard- Jones 
analogue in which the particles' coordinates vary continuously. A comparison 
between the two models allows us to discuss some exceptional, hardly realistic 
features of the original discrete system — which has been considered a prototype 
for nonequilibrium anisotropic phase transitions. 

3.1 Driven Lattice Gas 

The driven lattice gas [53) is a nonequilibrium extension of the Ising model with 
conserved dynamics. The DLG consists of a d-dimensional square lattice gas 
in which pair of particles interact via an attractive and short-range Ising-like 
Hamiltonian, 



Here Ck is the lattice occupation number at site k e Z'^ , and the sum runs over 
all the nearest-neig hbor (NN) sitefl Each lattice site has two possible states, 
namely, a particle (dk = 1) or a hole (cJk = 0) may occupy each site k. Cell 

^For d = 2, this corresponds with a lattice's connectivity 4. 




(3.1) 
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dimensions for the lattice are chosen so that the NN distance {bond length) was 
unity: |j— k| = 1 . A configuration is given by cr = {(Jk ; k G d-dimensional lattice}. 
We shall restrict ourselves to two-dimensional Ly x L_l lattices, that is, d = 2. 
Dynamics is induced by the competion between a heat bath at temperature 
T and an external driving field E which favors particle hops along one of the 
principal lattice directions, say horizontally (5i direction), as if the particles were 
positively charged. Consequently, for periodic boundary conditions, a nontrivial 
nonequilibrium steady state is set in asymptotically. This is formalized through 
a master equation similar to Eq. (j2.ip . where the transition probabilities per 
unit time w(cr cr') are the Metropolis ones, namely. 



Where E • 6 denotes the dot product between the field E (oriented along H direc- 
tion) and 5, which is the attempted particle displacement, i.e, 6 = (j — k)(ffj — 
cJk); and AH — H((j') — H((j). To be precise, w(ff — ^ a') stands for the rate 
for the particle-hole exchange (Jk ^ cjj when the configuration is cr; the par- 
ticle density p — N/(L||L_l) remains constant during the time evolution. The 
exchange ffk ^ corresponds to a particle jumping to a NN hole if Ck ^ crj; 
we assume w(cr — > cr') =0 otherwise, i.e., only NN particle- hole exchanges are 
allowed. This bias breaks detailed balance and establishes a nonequilibrium 
steady state. Typical configurations in the stationary regime of the DLG are 
shown in Fig. 12.11 

MC simulations by the biased Metropolis rate in Eq. p.2p reveal that, as 
in equilibrium, the DLG undergoes a second order phase transition. At high 
enough temperature, the system is in a disordered state while, below a critical 
point (at T < Te, where E = |E|) it orders displaying anisotropic phase segre- 
gation. That is, an striped rich-particle phase then coexists with its gas. It is 
also found that, for a lattice half filled of particles, the critical temperature Te 
monotonically increases with E from the Onsager valu^To — Tons — 2.269 Jkg^ 
to Too — L4Tons- This limit (E — > cxd) corresponds to a nonequilibrium critical 
point. As a matter of fact, it was numerically shown to belong to a universality 
class other than the Onsager one, e.g., MC data indicates Pdlg — 0-33 (instead 
of the Ising value (3 = 1/8 in two dimensions) for the order parameter critical 
exponent ^nj 151] [52] . Typical configurations for the DLG in the large field limit 
were shown in Fig. 12.11 

Statistical field theory is a complementary approach to the understanding 
of nonequilibrium ordering in the DLG. The derivation of a general mesoscopic 
description is still an open issue, however. Two different approaches have been 
proposed. The driven diffusive system (DDS) fl5 [l83[[51] . which is aphenomeno- 
logical Langevin-type equation aimed at capturing all the relevant symmetries, 
predicts that the current will induce a predominant mean-field behavior and, 
in particular, |3dds — V-2- The anisotropic diffusive system (ADS) f80], which 

■^See Appendix [Al 




(3.2) 
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Figure 3.1: Typical configurations during the stationary regime of the DLJF 
subject to a horizontal field of intensity E* = 1 . These graphs, which are for 
N = 900 particles and density p* — 0.30, illustrate, from left to right, (i) 
coexistence of a solid and its vapor (the configuration shown is for temperature 
T* — 0.20), (ii) Hquid-vapor coexistence (T* — 0.35), and (in) a disordered, 
fluid phase (T* — 0.50). The left-most graph shows a detail of the solid strip. 
The particles, which move in a square box of side -y/N/p*, are given here an 
arbitrary size, namely, diameter^ 1.1 cr. 



follows after a coarse graining of the master equation, rules out the relevance of 
the current and leads to the above-indicated MC critical exponent for E — > oo. 
However, the ADS approach reduces to the DDS for finite E, a fact which is 
hard to be fitted to MC data. In any case, field theoretical studies have con- 
stantly demanded further numerical efforts, and the DLG is nowadays the most 
thoroughly studied system showing an anisotropic nonequilibrium phase tran- 
sitions. The topic is not exhausted, however. On the contrary, there remain 
unresolved matters such as the above mentioned issues concerning critical and 
mesoscopic behaviors, and the fact that Te increases with E, which is coun- 
terintuitive |86[ I85j . Another significative question concerns the observation of 
triangular anisotropics at early times after a rapid MC quench from the ho- 
mogenous state (as if T = oo) to T< Te. The triangles happen to point against 
the field, which is contrary to the prediction from the DDS continuum equation 
ISIl- 

In this chapter we present a description of driven systems with continuous 
variation of the particles' spatial coordinates — instead of the discrete variations 
in the DLG — aiming for a new effort towards better understanding basic fea- 
tures of the DLG and its (and related) nonequilibrium phase transitions. This 
will provide a more realistic scheme for computer simulation of anisotropic fiu- 
ids. Our strategy to set up the model is to follow as closely as possible the 
DLG. That is, we analyze an off-lattice representation of the DLG, namely, a 
microscopically continuum with the same symmetries. Investigating these ques- 
tions happens to clarify the puzzling situation indicated above concerning the 
outstanding behavior of the DLG. 
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3.2 Driven Lennard-Jones fluid 

Consider a fluid consisting of N point-like particles in a two-dimensional L x L 
(L = L|| = Lj^) box with periodic (toroidal) boundary conditions. Interactions 
are according to a truncated and shifted Lennard-Jones (LJ) 12-6 potentialf] 

m-- 

V(r,j) = P^r(^^i^7^Lf^^^^' !^^^'<^'= (3.3) 

' \ 0, if Tij > Tc , 

where 

VLjW=4e[(cT/r)i2_(cT/r)'^] (3.4) 

is the full LJ potential. Here, Tij — \ri — rj | is the relative distance between 
particles i and j, e and ff are our energy and length units, respectively, and 
is the cut-off that we shall fix at — 2.5cr. The choice of this potential obeys 
to our strategy to set up the model following as close as possible the DLG (see 
Eq. p.ip V In fact, the potential defined in Eq. p. 31) is only one of the multiple 
possibilities for an attractive short-ranged potentiao The preferential hopping 
will be implemented as in the lattice, i.e., by adding a driving field to the poten- 
tial energy. Consequently, the familiar energy balance which enters in Eq. (13.21) 
is (assuming ke = 1 hereafter) H{c) = ^i<j where {continuum) config- 

urations c are here c = {ri,...,rN}. The particle attempted displacement is 
therefore 6 = r^ — rt. Lacking a lattice, the field is the only source of anisotropy, 
and any trial move should only be constrained by a maximum displacement in 
the radial direction. That is, we take < |5| < Smax) where Smax — 0.5cr in 
our simulations. We express the relevant quantities in terms of our energy e 
and length a scales, i.e., we introduce reduced unitf[£|. From these basic units, 
the temperature T, length L, number density p — N/L'^, internal energy U, and 
field E variables will be reduced respectively according to 

r = T/e , L* = L/(j , p* = pff2 , U* = U/e , and E* = Ea/e , (3.5) 

where * denotes reduced units. This model thus reduces for E to the well- 
known truncated and shifted LJ fluid, one of the most studied models in the 
computer simulation of fiuids [88l [89] . 



We studied this driven LJ fluid (DL JF) in the computer by the MC method 
using an extended "canonical ensemble" , namely, fixed values for N , p* , T* , and 
E*. Simulations involved up to N = 10^ particles with parameters ranging as 
follows: 0.5 < E* < 1.5, 0.20 < p* < 0.60, 0.15 < T* < 0.55. The typical config- 
urations one observes are illustrated in Fig. 13.11 As its equilibrium counterpart, 
the DLJF exhibits three different phases (at least): vapor, liquid, and solid (sort 



^ Often termed Weeks-Chandler-Andersen potential. 

*It has been shown that the LJ potential provides accurate results for a variety of fluids, 
e.g., Argon |88l . 

^This is very convenient in molecular simulations, and, most importantly, if we do not 
use reduced units, we might miss the equivalence between corresponding states {law of corre- 
sponding state). 
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Figure 3.2: Schematic pliase diagrams for the three models, as defined in the 
main text of this Chapter, showing ordered or striped (OP) and disordered (DP) 
phases. The left graph is for the DLG, for which To — Tons- The graph on the 
right is valid for both the NDLG, i.e., the DLG with next-nearest-neighbor 
(NNN) hops, for which To = 2.32 Tons, and the DLJF, for which Tq* = 0.459 



of close-packing phase; see the left-most graph in Fig. 13. ip . At intermediate 
densities and low enough T* , vapor and a condensed phase segregate from each 
other. The condensed droplet (see Fig. 13.11) is not near circular as it generally 
occurs in equilibrium, but strip-like extending along the field direction. A de- 
tailed study of each of these phases will be reported in Chapter 2] for a related 
system. 

A main observation is that the DLJF closely resembles the DLG in that both 
depict a particle current and the corresponding anisotropic interface. However, 
they differ in an essential feature, as illustrated by Fig. 13.21 That is, contrary 
to the DLG, for which the critical temperature Te increases with E, the DLJF 
shows a transition temperature Te which decreases with increasing E. The latter 
behavior was expectable. In fact, as E is increased, the effect of the potential 
energy in the balance Eq. p. 21) becomes weaker and, consequently, the cohesive 
forces between particles tend to become negligible. Therefore, unlike for the 
DLG, there is no phase transition for a large enough field, and Te — > for 
E — > CO in the DLJF. Confirming this, typical configurations in this case are 
fully homogeneous for any T under a sufficiently large field E. One may think 
of variations of the DLJF for which Te^oo — const > (see Chapter 2]), but 
the present one follows more closely the DLG microscopic strategy based on 
Eq. (123). 

Concerning the early process of kinetic ordering, one observes triangular 
anisotropies in the DLJF that point along the field direction. That is, the early- 
time anisotropics in the off-lattice case (right graph in Fig. l3.3p are similar to the 
ones observed in mesoscopic approaches (as, for instance, the DDS approach), 
and so they point along the field, contrary to the ones observed in the discrete 
DLG (left graph in Fig. [331). 
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Figure 3.3: Triangular anisotropics as observed at early times for E = 1 in com- 
puter simulations of the lattice (left) and the off-lattice (right) models defined 
in the main text. The field points along the 9. direction. The DLG configuration 
is for t = 6 X 10'* MC steps in a 128x128 lattice with N = 7372 particles and 
T = 0.4Tonsagor- The DLJF configuration is for t = 1 .5 x lO'^ MC steps, N = 10"^ 
particles, p* = 0.20 and T* = 0.23. 

3.3 Driven Lattice Gas with extended dynamics 

The above observations altogether suggest a unique exceptionality of the DLG 
behavior. This is to be associated with the fact that a driven particle is geo- 
metrically restrained in the DLG. In order to show this, we studied the lattice 
with an infinite drive extending the hopping to next -nearest-neighbors (NNN) 
[501 [nH [22] ■ That is, we extend interactions and accessible sites to the NNN. 
Thus, the sum in Eq. p.ip involves the eight next-nearest neighbors instead of 
the four nearest neighbors of the standard DLG. In this case, the particle-hole 
exchange dk Oj corresponds to a particle hop to a NNN hole if cr^ 7^ (Jj. 
As illustrated in Fig. 13. 4[ this introduces further relevant directions in the lat- 
tice, so that the resulting model, to be named here NDLG, is expected to behave 
closer to the DLJF. This is confirmed. For example, one observes in the discrete 
NDLG that, as in the continuum DLJF, decreases with increasing E — though 
from To = 2.32 Tons in this cas This is illustrated in Fig. 13.21 Specifically, 
we discuss this diflEerence between the DLG and the NDLG in Appendix [BJ in 
where we define an intermediate model which covers both situations. 

3.3.1 Correlations and Structure Factor 

There is also interesting information in the two-point correlation function and 
(equal-time) structure factor. The former, which measures the degree of order 
of the lattice, is defined for a half-filled lattice as 



where the steady average (■ • • ) involves averaging over j. Translational invariance 
is assumed. The structure factor, which is the Fourier transform of the two-point 

^For the Ising model with NNN interactions, one can easily derive theoretical estimates for 
the critical temperature I94| . 



C(i) = (cTjCTi+j)~1/4, 



(3.6) 
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Figure 3.4: Schematic comparison of the accessible sites a particle (at the center, 
marked with a dot) has for nearest-neighbors (NN) and next-nearest-neighbors 
(NNN) hops at equilibrium (left) and in the presence of a large horizontal field 
(right). The particle-hole exchange between neighbors may be forbidden (x), 
depend only on the potential energy (T), or occur with probability 1 (E). 




Figure 3.5: Surface plots of the two-point correlation (upper row) and the 
structure factor (lower row). These arc for the DLG (left column) and NDLG 
(right column) on a Ly x Lj^ = 128 x 128 lattice in the large field limit E = oo. 
Temperatures are T/Tons = 1 -58 and T/Tons = 1 -00 for the DLG and the NDLG, 
respectively. Notice the change of color code with each plot. 
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Figure 3.6: Parallel (squares) and transverse (triangles) components of the 
two-point correlation function above criticality for the DLG on a 128 x 128 
lattice with NN hops/interactions (filled symbols) at T = 1.58 and NNN 
hops/interactions (empty symbols) at T = 1. Temperatures normalized to Tons. 
The inset shows the x^^ power law decay in C|| for both discrete cases: DLG 
(o) and NDLG (x). 

correlation function, can be written in terms of the occupation variables as 




where k is a wave vector which has components K^^y = Inn-^^y/L with TLx,y = 
0, 1 , L — 1 . This is very useful to distinguish disordered configurations from 
inhomogeneous one. To avoid complications due to inhomogeneous ordered 
phases, we compute for both the DLG and NDLG two-body correlations and 
structure factors above criticality. We show in Fig. l3.5l the surface plots of C and 
its Fourier transform S for both dynamics. Analysis of the components (shown 
in Fig. 13. 6|) along the field, Cy = C(x, 0), and transverse to it, C_l = 0(0,^), 
shows that correlations are qualitatively similar for the DLG and the NDLG 
— although somewhat weaker along the field for NNN hops. That is, allowing 
for a particle to surpass its forward neighbor does not modify correlations. In 
fact, both the DLG and the NDLG display generic slow decay of two-point 
correlation^ [53]. It is generally accepted [TS] that the key ingredients which 
yield power-law decays are (i) a conserved dynamics, (ii) a non-equilibrium 
steady state, and (Hi) spatial anisotropy associated with the dynamics. The 
first ingredient alone (as in equilibrium) cannot lead to generic power laws due 
to the validity of the fluctuation dissipation theorem. The role of the second 

^This is in stark contrast with their equilibrium counterparts, the lattice gas (see Ap- 
pendix in which correlations decay exponentially (except at the critical point). The latter 
is a common feature of all the systems at thermal equilibrium. 
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Figure 3.7: Parallel (squares) and transverse (triangles) components of the struc- 
ture factor above criticality for the DLG on a 128 x 128 lattice for the DLG 
(filled symbols, left scale) at T = 1 .58 and for the NDLG (empty symbols, right 
scale) at T = 1 . Temperatures normalized to Tonsager- Here, nx,y — 128 K^^y /In. 



ingredient is to lift the theorem's constraint, so that the power laws are now 
generic. The role of the third one is more subtle, necessary only for producing 
power laws in the two-body correlations. 

The power-law behavioiH translates into a discontinuity of the structure fac- 
tor, namely, lim^^— soSy 7^ limK^j-jo Sj^, where Sy = S(Kx,0) and S± = S(0, Ky). 
This is clearly confirmed in Fig. 13.71 for both NN and NNN hops/interactions. 
Notice that this singularity, which do not appear under equilibrium conditions, 
is an indication of the dependence of the nonequilibrium steady state on the 
(anisotropic) dynamic. 



3.3.2 Critical behavior 



As outlined in Section 13.11 the derivation of a general mesoscopic description is 
still an issue under contention. In order to shed more light on the field theoretical 
description of driven fluids, we try to derive a Langevin-type equation for the 
NDLG. To this end, we firstly compute the order parameter critical exponent 
by means of standard finite size scaling techniques. Secondly, we derive the 
Langevin equation in order to describe the NDLG critical properties, by using 
the ADS approach [5D]. We also discuss the viability of other approaches, e.g., 
the DDS approach. 



Together with the Ising symmetry violation by E 0. 
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MC simulations 

Assuming a half-filled lattic(Jl we perform finite size scaling analysis for the 
NDLG by following the scheme proposed in [TQ consistent with the ADS theory 
[ST]. The order parameter is chosen as the structure factor 



ra = S(0,27t/U) , (3.8) 

(as suggested in [SSI ISZl ISH]), which carries the intrinsic anisotropics of the 
system. In order to perform systematic anisotropic finite size scaling we con- 
sidered system sizes 80 x 40, 25 x 50, 45 x 30, and 125 x 50. These aspect 
ratios satisfy fjj'"^^^" = 0.22 x Lj_, where Vj^/v|| « 1 /2 consistent with the ADS 
anisotropic spatial scaling [96l |97l |98] . A strong enough field is needed to avoid 
crossovers from the equilibrium regime. Notice also that, as we showed above 
for the NDLG, a saturating field suppresses the ordering, that is, Te = when 
E ^ oo. Therefore we choose an intermediate value for the field E = 12. The 
corresponding critical temperature is determined by using the Binder's fourth 
cumulant method [99], which is T12 — 1-45 ± 0.01. The obtained critical value 
T12 was employed for the finite size scaling analysis. 

Fig. 13. 81 show the scaled order parameter m upon the distance to the critical 
point A, for three different sizes of Ly. The former is rescaled by Lj'^^"'" whereas 

the latter is rescaled by The number of MC steps considered for each 

temperature was 160 millions. A perfect data collapse between the different 
system sizes is obtained by fixing vy w 1.25, y ~ ].22 and (3 ss 0.33, where vy, 
Y, and (3 are the corresponding critical exponents. Notice that these values are 
precisely the same for the DLG (with NN hops/interactions) [81]. As shown 
in Fig. 13.81 the asymptotic lines fit nicely with (3 sa 0.33 (upper line), and 
—y/2 « —0.61 (lower line). 

We have also computed the susceptibility, defined as the relative fluctuations 
of the order parameter. This reads 



X=-^((mV(-)^) (3.9) 

The rescaled susceptibility is shown in Fig. 13.91 The best collapse is obtained 
with the same values as before. Plotting the dimensionless Binder cumulant b 
as a function of the rescaled distance to the critical point with vy w 1 .25, again, 
nearly perfect collapse is obtained for all system sizes (as depicted in Fig. l3.10() . 

This set of MC results support that both the DLG and the NDLG belong to 
the same universality class, i.e., both dynamics yield the same critical indexes. 
In particular, it follows that the critical properties remains unchanged when 
extending hops and interactions from the nearest-neighbors to the next -nearest- 
neighbors. 



^By following simple symmetry arguments, the eritieal density is N/LxLii = 1/2. 



3.3 Driven Lattice Gas with extended dynamics 



29 




-2-10 12 3 



Figure 3.8: Log-log plot of the order parameter rescaled by Lj*^^^" vs. ALy^^", 
for different sizes: (+) 125 x 50, (*) 45 x 30, and (o) 80 x 40. The upper straight 
line is for P = 1 /3, whereas the lower one indicates y — 0.3. 



5.5 
4.5 



> 3.5 



2.5 



0.5 







1 


















































J 








I 
























* 












o' 












1 






1 





XL 



l/v„ 



Figure 3.9: Log-log plot of the susceptibility rescaled by L..^^^^^ vs. AL|/^" 



Symbols are as in Fig. [37 
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Figure 3.10: Scaling plot of the fourth Binder cumulant b vs. /M-^i^^" . Symbols 
are as in Fig. 



Mesoscopic description 

The ADS approach lias been employed successfully to describe the critical be- 
havior of the DLG [8CT|. The resulting Langevin-type equation, derived from 
a coarse-graining process from the DLG master (microscopic) equation sets up 
the anisotropy as the main noncquilibrium ingredient and leads to the right 
critical exponent in the large field limit (|3 « 1/3). The associated Langevin 
equation for the coarse-grained density 4)(r, t) reads 

9t4)(r,t) = -Vl4) + (T+ a)Vi4) + ^^l<t>' + f VfcJ) + £,(r,t) (3.10) 

Here, the last term £, stands for a Gaussian white noise, i.e., (f,(r, t)) = and 
(f,(r, t)£,(r', t')) = 6(r — r')5(t — t'), representing the fast degrees of freedom, 
and T, a and g are model parameters. 

This approach is certainly valuable because it represents a detailed connec- 
tion between microscopic dynamics and their mesoscopic descriptions. There- 
fore, this important feature suggests that its validity can be extended to describe 
the critical behavior of the NDLG, which was studied by computational means 
in the previous subsection. Such a description include the microscopic details of 
diagonal dynamics which, in particular, should allow us to distinguish between 
the DLG and the NDLG. This is in contrast of other more heuristic approaches, 
as the DDS [TS] [HSl HI]: which aim for the main symmetries of the DLG. The 
DDS approach is the natural extension of the conserved cf)^ theory for the Ising 
equilibrium model |100j . The most relevant prediction of the DDS is the mean 
field behavior (3 = 1 /2, which conflicts with MC simulations. This is not sur- 
prising, because the DLG behavior is more complicated than its equilibrium 
counterpart, in which the symmetries and the influence of the dynamics is well 
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established. The same situation occurs when extending dynamics to the NNN: 
the DDS (heuristic) scheme leads again to the mean-field order parameter crit- 
ical exponent. As a matter of fact, the DDS approach can not account for 
the diagonal (microscopic) degrees of freedom. Therefore, the DDS Langevin 
equation does not describe the NDLG critical properties. 

The question is whether the ADS approach describes properly the critical 
behavior observed in the NDLG. This is addressed in this section. 

Next, we outline the derivation of a Langevin equation for the NDLG by 
following the scheme proposed in |80) . (this derivation can be found in greater 
detail in Appendix [C]) . 

Let us define a density variable 4)(r, t) averaged in a region of volume v over 
which the original microscopic occupation variable (cTk in Eq. ([33J) were aver- 
aged. The system evolves from a (coarse-grained) configuration y to another 
y' by choosing randomly a particle at point r and exchanging it with one of 
its next nearest neighbors in the direction \i [80]. Notice that, in contrast with 
the previous ADS approach for the DLG with NN interactions, now |j. involves 
the known parallel || and transversal _L directions as well as the two additional 
directions which allow for the diagonal jumps / and \. If v is large enough, 
then 4) can be considered a continuum function of r so we have 

y' = {4) +v"^ V^5(r' -r), where cj) G y} (3.11) 

The distribution P(y, t) accounts for the statistical weight of each configuration 
y at time t and evolves according to the following Markovian master equation 

DP 



dr my ^ y )P(y, t] - a(y' ^ y)P(y', t) 



9tP{y,t)= Y. 'i^f(^) 

n={ii,_L,/,\r 

(3.12) 

where 0.(y — > y') stands for the transition probability per unit time (transition 
rate) from y to y' and f (ri) is an even function which account for the amount of 
mass attempted to be displaced. That is, Eq. (|3.11l) - (|3.12p represent a proccess 
in which a configuration y is exchanged with the infinitesimal neighbor of r (y') 
in the |j. direction. As usual, the transition rates are taken to be a product of 
a function of the entropy times the same function of the difference between the 
configurations plus a term due to the effect of the driving field [351 , namely, 

a[y ^ y') = D(AS{ti)) • D{AH{ti) + H,). (3.13) 

Here, = r[\L - ^[^ — cf)^) + 0[v^^), and D is a function satisfying detailed 
balance. Importantly, in the Langevin-type equations framework, e represents 
a coarse-grained field acting upon a given configuration rj. Indeed, the detailed 
dependence of the mesoscopic coefficients on the microscopic field introduced 
in Eq. (13. 2p remain still as an open issue [llj . This constraint ensures that 
the limiting case of e = the steady state solution of the master equation is 
the canonical one, i.e. P(y,t = 00)1^=0 = e-"t*l. At a mesoscopic level, the 
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equivalent to Eq. (|3.1|) is the standard c})^ (Ginzburg-Landau) Hamiltonian, 
and following a similar procedure as here the equilibrium Langevin equation is 
the so called model B |100j . The structure of the free energy consist of two 
contributions: a entropic and a energetic one which are given by, respectively 



S(Tl) 

H(ti) = 



dr 



4!^ 



dr ( 1 [V^f + I 



(3.14) 



Here a and g are entropic coefficients whereas the parameter t comes from the 
energetic contribution. From Eq. (|3.12l) a Fokker-Planck equation is derived 
by means of a Kramers-Moyal fT] expansion. The derivation is similar to the 
one in Ref. [80], but taking into account the transversal degrees of freedom 
for the NDLG. After using standard techniques in the theory of stochastic pro- 
cesses [T] we derive from the Fokker-Planck equation its stochastically equivalent 
Langevin equation using the Ito prescription, which reads 



where 

Af, = H.£(1 -4)2), 



qiZl,Z2j 



(3.15) 



drrif(ri)D(rizi)D(TiZ2), 



drTi2f(Ti)D{Tizi)D(riZ2) 



(3.16) 



The time has been rescaled by v^^ and f,|x(r, t) is a delta-correlated Gaussian 
white noise in the \i direction. After some algebra, the terms depending on both 
diagonal components / and \ can be split into the two parallel and transversal 
to the field components. We now focus on the critical region and discard irrel- 
evant terms in the renormalization group sense by naive power counting. We 
perform the following anisotropic scale transformations: 



t T^^t , 



(3.17) 



4> 



Expanding the Langevin equation Eq. p.lSp around t = keeping only the 
more relevant terms, we found that the leading terms of the theory are V^(|) and 
V?4) which, together with imposing invariance of the time scale, the transverse 
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spatial interaction, and the transverse noise, under scale transformations, lead 
to s =2, 6 = (s + d — 3]/2, and z = 4. According to this and after some 
cumbersome algebra, the Langevin equation for finite (coarse-grained) driving 
field contains a large bunch of field-dependent terms. Similarly to what occurs 
in the ADS approach for the DLG, just taking a large enough value of the field 
most of those terms vanish. By doing this, one is led to the following Langevin 
equation 

9t4)(r, t) = -Vl4) + (t + ^a)\/i4> + + aVf cf) + t(r, t) (3.18) 

The resulting equation resemble the one for NN interactions (see Eq. I|3.10p ). 
although in Eq. (|3.18l) appears additional prefactors of entropic nature. This is 
consistent with the fact that the model with the extended dynamics, i.e., the 
NDLG, possesses more neighbors than the DLG. The present equation for NNN 
interactions contains all the same symmetries as Eq. p.lOp and both lead to the 
same critical behavior. The value of the order parameter critical exponent that 
follows from a renormalization group analysis of Eq. (|3.18p is (3 = 1/3, which is 
consistent with the observed in the MC simulations described above. 

3.4 Conclusions 

The above observations show that the nature of correlations is not enough to 
determine the phase diagram. There are already indications of this from the 
study of equilibrium systems, and also from other nonequilibrium models. That 
is, one may have a qualitatively different phase diagram but essentially the 
same two-point correlations by modifying the microscopic dynamics. It follows, 
in particular, that the exceptional behavior of the DLG cannot be understood 
just by invoking the two-point correlation function C(i) — or, equivalently, the 
structure factor S(k) — or crude arguments concerning symmetries. The fact 
that particles are constrained to travel precisely along the two principal lattice 
directions in the DLG is the cause for its singular behavior. Or in other words, 
the lattice geometry acts more efficiently in the DLG as an ordering agent than 
the field itself, which occurs rarely in actual cooperative transport in fluids. 
This means that the ordering in the DLG is a geometric effect rather than 
dynamic. Allowing jumping along intermediate directions, as in both the NDLG 
and the DLJF, modifies essentially the phase diagram but not features — such 
as power-law correlations — that seem intrinsic of the nonequilibrium nature 
of the phenomenon. This special arrangement of particles is also lacking, for 
instance, when other orientations of the driving field |101j are considered. This 
also would explain certain irregular configurations observed in three dimensions 
[TT| . since this constraint is also lifted in the d = 3 case. 

Rutenberg and Yeung [92j also performed quenching experiments for varia- 
tions of the DLG0. They showed, in particular, that minor modifications in the 



They allow for NNN hops, but the Hamiltonian only involves the NN interactions. 
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DLG dynamics may lead to an inversion of the triangular anisotropies during 
the formation of clusters which finally condense into strips. Our observations 
above prove that such nonuniversal behavior goes beyond kinetics, namely, it 
also applies to the stationary state. 

The unique exceptionality of the DLG has some important consequences. 
One is that this model involves features that are not frequent in nature. There 
are situations in which a drive induces stripes but not necessarily DLG behavior 
[ST]. The fact that a particle impedes the freedom of the one behind to move 
along a for a large enough field may only occur very seldom in cooperative trans- 
port in fluids. The NDLG is more realistic in this sense. In any case, the great 
effort devoted to the DLG during two decades has revealed important proper- 
ties of both nonequilibrium steady states and anisotropic phase transitions. On 
the other hand, it ensues that, due to its uniqueness, the DLG does not have 
a simple off-lattice analog. This is because the ordering agent in the DLG is 
more the lattice geometry than the field itself The fact that one needs to be 
very careful when modeling nonequilibrium phenomena — one may induce both 
a wrong critical behavior and an spurious phase diagram — ensues again in this 
example. This seems not to be so dramatic in equilibrium where, for example, 
the lattice gas is a useful oversimplification of a LJ fiuid. 

Concerning criticality, according to our observations we found that the aniso- 
tropic diffusive system approach include the microscopic details of transverse 
dynamics which, in particular, allow one to distinguish between the DLG and 
the NDLG. This is in contrast with other, more heuristic approaches, e.g., the 
driven diffusive system. The ADS approach lead to the right order parameter 
critical exponent given by MC simulations. We analyzed the order parameter, 
magnetization, and the fourth Binder cumulant in order to compute the critical 
exponents. By using anisotropic finite size scaling we have shown that, in spite 
of structural differences between them, both models belong to the same uni- 
versality class, which is characterized by an order parameter critical exponent 
P«l/3. 

Consequently, their critical properties do not depend on this particular ex- 
tension of the dynamics. There are other nonequilibrium models in which the 
corresponding critical behavior changes by such an extension of the dynamic 
from NN to NNN. That is, for instance, the case of a closely related lattice sys- 
tem (with NN exclusion under a driving field) devised by Dickman [93j : When 
one extends the dynamics to NNN [91] , the nature of the order-disorder transi- 
tion which takes place in this model changes dramatically. We shall discuss on 
this model further in the next chapter. 

Finally, we remark that in the novel fluid model introduced in this chapter, 
the particle 'infinite freedom' is realistic, as it is also the LJ potential. It may 
contain some of the essential physics in a class of nonequilibrium anisotropic phe- 
nomena and phase transitions. On the other hand, the model is simple enough 
to be useful in computer simulations, and it is endowed of even-simpler and 
functional lattice analogs such as the NDLG. In the following chapter (Chap- 
ter |4]) we shall develop in greater detail this matter. 



Chapter 4 



Prototypical model for 
anisotropic phenomena in 
nonequilibrium fluids 

In this chapter, we describe short-time kinetic and steady-state properties of 
the non-equihbrium phases, namely, (mono- and polycrystalhne) sohd, hquid 
and gas anisotropic phases in a novel driven Lennard- Jones fluid. This is a 
computationally-convenient model which exhibits a net current and striped 
structures at low temperature, thus resembling many situations in nature. We 
here focus on both critical behavior and some details of the early nucleation 
process, by far less studied than later regimes. In spite of the anisotropy of the 
late-time "spinodal decomposition" process, earlier nucleation seems to proceed 
by Smoluchowski coagulation and Ostwald ripening, which are known to account 
for nucleation in equilibrium, isotropic lattice systems and actual fluids. On the 
other hand, a detailed analysis of the system critical behavior rises some in- 
triguing questions on the role of symmetries; this concerns the computer and 
field-theoretical modeling of nonequilibrium fluids. 

4.1 Motivation 

The microscopic mechanism of stripe formation in particular materials is still 
not cleai0. Hence, lacking theory for the "thermodynamic" instabilities causing 
the observed striped structures, one tries to link them to the microscopic dy- 
namics of suitable model systems. As we explained in the previous chapters, the 
driven lattice gas (DLG) |53) . namely, a model system in which particles diffuse 
under an external driving field, has been up to now a theoretical prototype of 
anisotropic behavior [TTJ [TSl [2D] ■ This model was shown in Chapter [3] to be un- 



^For instance, in high— temperature superconductors 11021 . questions such as whether or 
not the Coulomb interactions are important in these kind of materials remain unresolved. 
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Figure 4.1: The place of our novel off-lattice models (as, for instance, the DLJF 
introduced in the previous chapter) in the realm of computer simulation of fluids. 



realistic in some essential sense, however. Particle moves in the DLG are along 
the principal lattice directions, and any site can hold one particle at most, so 
that a particle impedes the one behind to jump freely along the direction which 
is favored to model the action of the field. Consequently, the lattice geome- 
try acts more efficiently in the DLG as an ordering agent than the field itself, 
which occurs rarely — never so dramatically — in actual cooperative transport 
in fluids. In fact, actual situations may in principle be more closely modeled 
by means of continuum models, and this peculiarity of the DLG implies that it 
lacks a natural off-lattice extension. 

Here we present and analyze numerically a novel nonequilibrium off-lattice, 
Lennard- Jones (LJ) system which is a candidate to portray some of the anisotro- 
pic behavior in natur(|^. The model, which involves a driving field of intensity 
E, reduces to the celebrated (equilibrium) LJ fluid [88l |89] as E — ) 0. For any 
E > 0, however, it exhibits currents and anisotropic phases as in many ob- 
servations out of equilibrium. In particular, as the DLG, our model in two 
dimensions shows striped steady states below a critical point. We also observe 
critical behavior consistent with the equilibrium universality class. This is rather 
unexpected in view of the criticality reported both for the DLG. On the other 
hand, concerning the early-time relaxation before well-defined stripes form by 
spinodal decomposition, we first observe — as in previous studies of relaxation 
towards equilibrium — effective diffusion of small droplets, which is followed by 
monatomic diffusion probably competing with more complex processes. It is 
very likely that our observations here concerning nucleation, coexistence, criti- 
cality, and phases morphology hold also in a number of actual systems. 

Before we proceed, it is worth mentioning that the off-lattice model^ devised 
over the course of this thesis fills a void in the physics of nonequilibrium fiuids. 

^This is in fact a related model to the driven Lennard- Jones fluid (DLJF) introduced in 
the previous Chapter. 

''Here we include the DLJF model defined in Chapter |3] 
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This is schematized in Fig. 14.11 As we explained in Chapter [1] lattice models 
are convenient, simplified models which often capture the main ingredients of 
equilibrium systems, particularly, close to the critical points. For instance, the 
lattice gas may be a useful representation of, e.g., a Lennard- Jones fluid. How- 
ever, rather surprisingly, nonequilibrium critical phenomena are by far more 
widely studied in nonequilibrium extensions of lattice models than in their con- 
tinuum counterparts. In this chapter, we provide a more realistic — without loss 
the efficiency of lattice models — , convenient model for the study of nonequi- 
librium phase transitions in fluids, and we suggest simple procedures to extend 
the study of actual fluids to the nonequilibrium realm. 

In section l42l we define the off-lattice model, and section l473l is devoted to 
the main results as follows. § I4.3.1l describes the early-time segregation process 
as monitored by the excess energy, which measures the droplets surface. § 14.3.21 
describes some structural properties of the steady state, namely, the radial and 
azimuthal distribution functions, and the degree of anisotropy. § 14.3.31 depicts 
some transport properties. And § 14. 3. 41 is devoted to an accurate estimate of the 
liquid-vapor coexistence curve and the associated critical parameters. Finally, 
Section [4.41 contains the main conclusions and final comments. 



4.2 The model 

Consider N particles of equal mass (set henceforth to unity) in a d— dimensional 
square box, L'^, with periodic boundary conditions. Interactions are via the 
truncated and shifted pair potential [88] (already defined in Chapter [3] see 
Eq. (1331)): 

*i^> i;;;::, 

where 

4)Lj(r)=4e[((T/r)i2_(^/^^6] (42) 

is the full LJ potential, r stands for the interparticle distance, and tc is a cut- 
off that we set — 2.5cr. The parameters cr and e are, respectively, the 
characteristic length and energy — that we use in the following to reduce units 
as in Eq. ([33]) . 

Time evolution is by microscopic dynamics according to the transition prob- 
ability per unit time (rate): 

a)(t' (Tl^n') ^x'"" xmin{l,e-^*/T}. (4.3) 

Here, x'^' j which is a nonequilibrium extension of the Metropolis rate |5T, reads 



X'^' = ^[1 + tanh(E«-5)], 



(4.4) 
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Figure 4.2: Schematic representation of the region (grey) which is accessible to 
a given particle as a consequence of a trial move for E = (left-handside) and 
E = 00 (right-handside), assuming the "infinite" field points % horizontally. 

E is the intensity of a uniform external field along the horizontal direction, say 
■^5 T) = {ri , . . . , Tim} stands for any configuration of energy 



where rt is the position of particle i that can be anywhere in the d— torus, r)' 
equals ri except for the displacement of a single particle by 6 = r- — rt, and 
AO = (t>[r\'] — 0(ri) is the cost of such displacement. 

It is to be remarked that x'^') as defined in (|4.4p . contains a drive bias 
(see Fig. 14. 2p such that the rate (|4.3|) lacks invariance under the elementary 
transitions ri ^ r\\ Consequently, unlike in equilibrium, there is no detailed 
balance for toroidal boundary conditions if E > 0. 

Here we describe the results from a series of Monte Carlo (MC) simulations. 
In order to reduce unnecessary interparticle distance calculations, we used a 
standard neighbor-list algorithm [88]. Our simulations concern fixed values of 
N, with N < 10'*, particle density p = N/L'^ within the range p e [0.2, 0.6] , and 
temperature T e [10^^,10^] . Following the fact that most studies of striped 
structures, e.g., many of the ones are mentioned in the Chapter [2l concern two 
dimensions — in particular, the DLG critical behavior is only known with some 
confidence for d — 2 [TTJ [81] [82] — we restricted ourselves to a two dimensional 
torus. We fixed the maximum particle displacement is 6max — 0.5 in our simu- 
lations as shown in Fig. 14.21 We deal below with steady-state averages over 10^ 
configurations, and kinetic or time averages over 40 or more independent runs. 

The distribution of displacements 6 is uniform, except that the new particle 
position is (most often in our simulations) sampled only from within the half- 
forward semi-circle of radius S^ax centered at r,,, as illustrated on the right-hand 
graph of Fig. 1121 This is because the infinite-field limit, E — > oo, turns out to be 
most relevant, and this means, in practice, that any displacement contrary to the 
field is forbidden. This choice eliminates from the analysis one parameter and, 
more importantly, it happens to match a physically relevant case. As a matter 
of fact, simulations reveal that any external field E > induces a flux of particles 
along 9, — which crosses the system with toroidal boundary conditions — that 
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monotonically increases with E, and eventually saturates to a maximum. This 
is a realistic stationary condition in which the thermal bath absorbs the excess 
of energy dissipated by the drive. Our simulations here concern only E = oo in 
order to maximize the nonequilibrium effecttlfl. 

4.3 Main results 

Figure illustrates late-time configurations, i.e., the ones that typically char- 
acterize the steady state, as the temperature T is varied. These graphs already 
suggest that the system undergoes an order-disorder phase transition at some 
temperature Te. This happens to be of second order for any E > 0, as in the 
equilibrium case E = 0. We also observe that decreases monotonically with 
increasing E, and that it reaches a well-defined minimum. Too, as E — > oo. 

Figure also shows that, at low enough temperature, an anisotropic in- 
terface forms between the condensed phase and its vapor; this extends along 9. 
throughout the system at intermediate densities. 

4.3.1 Phase segregation kinetics 

Skipping microscopic details, the kinetics of phase segregation at late times 
looks qualitatively similar to the one in other nonequilibrium cases, including 
the DLG 103 and both molecular-dynamic [104; and Cahn-Hilliard |105j rep- 
resentations of sheared fiuids, while it essentially differs from the one in the 
corresponding equilibrium system. This is illustrated in Fig. 14.41 One observes, 
in particular, condensation of many stripes — as in the graph for t = 10^ in 
Fig. 14.41 — into a single one — as in the first three graphs at the bottom row 
in Fig. 14.31 This process corresponds to an anisotropic version of the so-called 
spinodal decomposition 106 , which is mainly characterized by a tendency to- 
wards minimizing the interface surface as well as by the existence of a unique 
relevant length, e.g., the stripe width |103j . We focus here on a detailed study of 
the early-time separation process instead of the late regime, which has already 
been studied for both equilibrium |107[|108] and nonequilibrium cases, including 
the DLG [1031 1109] . This is because detailed descriptions of early nonequilib- 
rium nucleation are rare as compared to studies of the segregation process near 
completion. 

Following an instantaneous quench from a disordered state into T< Too owe 
observes in our case that small clusters form, and then some grow at the expenses 
of the smaller ones but rather independently of the growth of other clusters of 
comparable size. This corresponds to times t < 10^ in Fig. 14. 4[ i.e., before 
many well-defined stripes form. We monitored in this regime the excess energy 
or enthalpy, H (t) , measured as the difference between the averaged internal 
energy at time t > and its stationary value. This reflects more accurately the 
growth of the condensed droplets than its size or radius, which are difficult to 

^We are interested in the far-from-equilibrium effects induced by the external driving field 
rather than in weak deviations from equilibrium. 
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and p — 0.30. The field is oriented along the horizontal direction. 
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Figure 4.4: Typical configurations for E = (top row) and E — > oo (bottom row) 
as time proceeds during relaxation from a disordered state (as for T ^ oo) at 
t = 0. The time here is proportional to the MC steps (MCS) as indicated. This 
is for N = 7500, p = 0.35, and T = 0.275, below the corresponding transition 
temperature. 
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Figure 4.5: Time evolution of the enthalpy per particle for N — 7500, p = 0.35 
and, from top to bottom, T = 0.200, 0.225, 0.250, and 0.275. Straight lines are 
a guide to the eye; the slope of each line is indicated. The inset shows the detail 
at early times. (For clarity of presentation, the main graph includes a rescaling 
of the time corresponding to the data for T — 0.250, 0.225 and 0.200 by factors 
2, 3 and 3, respectively.) 

be estimated during the early stages |110[ llllj . In addition to its theoretical 
importance, the excess of energy may be relevant for experimentalists. This is 
because H (t) may be determined in microcalorimetric experiments [112] . 

The time development of the enthalpy density h. (t] = H (t] /N is depicted 
in Fig. 14.51 This reveals some well-defined regimes at early times. 

The first regime, (a) in the inset of Fig. l475l follows a power law t^^ with « 
1 / 6 — which corresponds to the line shown in the graph — independently of the 
temperature investigated. This is the behavior predicted by the Smoluchowski 
coagulation or effective cluster diffusion [99j . The same behavior was observed in 
computer simulations for E = and also reported to hold in actual experiments 
on binary mixtures |110[ 1112) . Surprisingly, this suggests the early dominance 
of a rather stochastic mechanism, in which the small clusters rapidly nucleate, 
which is practically independent of the field, i.e., it is not afi'ected in practice by 
the drive. It is also noticeable that, in contrast with simulations on lattice gases 
in which this regime is rather noisy, a slope t^^/^ appears clearly in Fig. 14.51 
The indication of some temperature dependence in equilibrium [111] , which is 
not evident here, might correspond to the distinction between deep and shallow 
quenches made in Ref. jllO) that we have not investigated out of equilibrium. 

At later times, there is a second regime, (b) in Fig. 14. 5[ in which the 
anisotropic clusters merge into filaments and, finally, stripes. We observe in this 
regime that 6 varies between 0.3 and 0.6 with increasing T. Ostwald ripening 
|113) , consisting of monomers difi^usion, predicts 9 = 1/3 for fluids at thermal 
equilibrium. The situation here in our nonequilibrium is more complicate: it is 
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Figure 4.6: Details of the structure in the low— T, solid phase as obtained by 
zooming into configurations such as the ones in Fig. 14.31 This is for T — 0.05, 
0.12, and 0.25, from left to right, respectively. 



likely that regime (b) describes a crossover from a situation which is dominated 
by monomers at low enough temperature to the emergence of other mechanisms 
|108| 1114] which might be competing as T is increased. 

Finally, one observes a regime, (c) in Fig. 14.51 which corresponds to the 
beginning of spinodal decomposition. 

4.3.2 Structure of the steady state 

Concerning the local structure in the stationary regime, for any E > 0, the 
anisotropic condensate changes from a solidlike hexagonal packing of particles 
at low temperature (e.g., T = 0.10 in Fig. 14.31) . to a polycrystalline or perhaps 
glass-like structure with domains which show a varied morphology at, e.g., 
T = 0.12. The latter phase further transforms, with increasing temperature, 
into a fluid-like structure at, e.g., T — 0.30 and, finally, into a disordered, 
gaseous state. 

More specifically, the typical situation we observe at low temperature is 
illustrated in Fig. 14.61 At sufficiently low temperature, T = 0.05 in the example, 
the whole condensed phase orders according to a perfect hexagon with one of 
its main directions along the field direction 9.. This is observed in approximately 
90% of the configurations that we generated at T = 0.05, while all the hexagon 
axes are slanted with respect to it in the other 10% cases. As the system is 
heated up, the stripe looks still solid at T = 0.12 but, as illustrated by the 
second graph in Fig. 14. 6[ one observes in this case several coexisting hexagonal 
domains with different orientations. The separation between domains is by 
line defects and/or vacancies. Interesting enough, as it will be shown later on, 
both the system energy and the particle current are practically independent of 
temperature up to, say T = 0.12. The hexagonal ordering finally disappears in 
the third graph of Fig. 14.61 which is for T — 0.25; this case corresponds to a 
fluid phase according to the criterion below. 

A close look to the structure is provided by the radial distribution (RD), 
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Figure 4.7: Data from simulations for N = 7000 and p = 0.35. The main grapli 
shows the degree of anisotropy, as defined in the main text, versus temperature. 
The vertical dotted line denotes the transition temperature. The lower (upper) 
inset shows the radial (azimuthal) distribution at T = 0.20, full line, and T — 

0. 30, dashed line. 

1. e., the probability of finding a pair of particles a distance r apart, relative to 
the case of a random spatial distribution at same density. This is shown in 
the lower inset of Fig. 14.71 At fixed T, the driven fluid is less structured than 
its equilibrium counterpart, suggesting that the field favors disorder. This is 
already evident in Fig. 14.31 and it also follows from the fact that the critical 
temperature decreases with increasing E. 

Further information on the essential anisotropy of the problem is beared by 
the azimuthal distribution (AD), which accounts for the longe range ordering. 
This is defined 



where 0ij G [0, Ztt) is the angle between the line connecting particles i and j 
and the field direction Except at equilibrium, where this is uniform, the AD 
is 7t/2— periodic with maxima at Ictt and minima at 1c7t/2, where k is an integer. 
This means that the field favors long-range longitudinal ordering instead of 
order along other directions. This is consistent with the observations in which 
the orientations of the hexagons in the solidlike phase are parallel to the field 
direction. The AD is depicted in the upper inset of Fig. 14.71 

We also monitored the degree of anisotropy, defined as the distance D 




(4.7) 



D = 



de |a(e)-i|, 



(4.8) 







which measures the deviation from the equilibrium, isotropic case, for which 
a(6) — 1, independent of 9. The function (|4.8p . which is depicted in the main 
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graph of Fig. 14.71 reveals the existence of anisotropy even above the transition 
temperature. This shows the persistence of nontrivial two-point correlations 
at high temperatures which has been demonstrated for other nonequilibrium 
models [55] . 

4.3.3 Transport properties 

The current is highly sensitive to the anisotropy. The most relevant information 
is carried by the transverse-to-the-field current profile jj^, which shows the dif- 
ferences between the two coexisting phases (Fig. 14. 8p . Above criticality, where 
the system is homogeneous, the current profile is flat on the average. Otherwise, 
the condensed phase shows up a higher current (lower mean velocity) than its 
mirror phase, which shows up a lower current (higher mean velocity). Both the 
transversal current and velocity profiles are shown in Fig. 14.81 The current and 
the density vary in a strongly correlated manner: the high current phase cor- 
responds to the condensed (high density) phase, whereas the low current phase 
corresponds to the vapor (low density) phase. This is expectable due to the 
fact that there are many carriers in the condensed phase which allow for higher 
current than in the vapor phase. However, the mobility of the carriers is much 
larger in the vapor phase. The maximal current occurs in the interface, where 
there is still a considerable amount of carriers but they are less bounded than 
in the particles well inside the bulk and, therefore, the field drives easily those 
particles. This enhanced current effect along the interface is more prominent in 
driven lattice models (notice the large peak in the current profile in Fig. 14. 8p . 
Thus, these are clear cases of interfacially driven transport |115| . Moreover in 
both lattice cases, DLG and NDLG (already studied in Chapter there is 
no difference between the current displayed by the coexisting phases because 
of the particle-hole symmetry. Such a symmetry is derived from the Ising-like 
Hamiltonian in Eq. p.ip and it is absent in the off-lattice case. 

Regarding the comparison between off-lattice and lattice transport proper- 
ties. Fig. 14.91 shows the net current j, defined as the mean displacement per 
MC step per particle, as a function of temperature. Saturation is only reached 
at jraax = 46Taax/37t whcu T — > oo. The current approaches its maximal 
value logarithmically, i.e., slower than the exponential behavior predicted by 
the Arrhenius law. The way this limit is approached is illustrated in the inset of 
Fig. 14. 91 The sudden rising of the current as T is increased can be interpreted as a 
transition from a poor-conductor (low-temperature) phase to a rich-conductor 
(high-temperature) phase, which is reminiscent of ionic currents This 
behavior of the current also occurs in the DLG. Revealing the persistence of 
correlations, the current is nonzero for any, even low T, though it is small, and 
roughly independent of T, in the solid phase. 

An insight of the transition points between the different phases is obtained 
from the temperature dependence of the mean potential energy per particle 
f = N^^ (3'('n)) a-nd from the net current j. In particular, one may estimate 
the transition point as the condensed strip changes from solid to liquid (T « 
0.15) and finally changes to a fully disordered state (T ~ 0.31), i.e., disordered 
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Figure 4.8: . Transverse-to-the-field current profiles below criticality. The 
shaded (full) line corresponds to the current (velocity) profile of the off-lattice 
model. For comparison we also show the current profile of the DLG (circle- 
dotted line). Since each distribution is symmetric with respect to the system 
center of mass (located here at L/2) we only show their right half parts. 

state. This is shown in Fig. 14.91 which exhibits well-defined changes of slope in 
both magnitudes when the phases transform. Notice also that the mean energy 
behaves linearly with temperature for T e (0.12,0.3) in the liquid phase. 

4.3.4 Coexistence curve 

One of the main issues concerning the steady state is the liquid-vapor coexis- 
tence curve and the associated critical behavior. The (nonequilibrium) coexis- 
tence curve may be determined from the density profile transverse to the field. 
This is illustrated in Fig. QUI 

At high enough temperature — in fact, already at T = 0.35 in this case for 
which the transition temperature is slightly above 0.3 — the local density is 
roughly constant around the mean system density, p = 0.35 in Fig. 14.101 As 
T is lowered, the profile accurately describes the existence of a single stripe 
of condensed phase of density p+ which coexists with its vapor of density p_ 
(p_ < p+). The interface becomes thinner and smother, and p+ increases while 
p_ decreases, as T is decreased. 

As in equilibrium, one may use p+ — p_ as an order parameter. The result 
of plotting p+ and p_ at each temperature results in the non-symmetric liquid- 
vapor coexistence curve shown in Fig. 14.111 The same result follows from the 
current, which in fact varies strongly correlated with the local density. Notice 
that the accuracy of our estimate of p± is favored by the existence of a linear 
interface, which enable us to get closer to the critical point than in equilibrium. 
Lacking a thermodynamic theory for "phase transitions" in nonequilibrium liq- 
uids, other approaches have to be considered in order to estimate the critical 
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Figure 4.9: Temperature dcipendenee of the mean energy (squares; the scale is on 
the right axis) and normalized net current (circles; scale at the left) for N = 5000 
and p = 0.30 under "infinite" field. The inset shows the T— dependence of the 
current over a wider range. 




Figure 4.10: Density profiles transverse to the field for N = 7000, p = 0.35, and 
different temperature, as indicated. The coexisting densities, p±, are indicated. 



48 



Anisotropic phenomena in nonequilibrium fluids 



0.35r 



0.3- 



a25h ° 

a 



-1 ■ 1 — 

fluid phase 



liquid - •uiipor 



a 



"0" 



0.2 



0.4 
P 



0.0 



0.8 



Figure 4.11: Coexistence curve (squares) for the LJ nonequilibrium model ob- 
tained from the density profiles in Fig. 14.101 The fluid phase and the coexistence 
region are indicated. The triangles are the arithmetic mean points, which serve 
to compute the critical parameters. The large circle at the top of the curve 
locates the critical point, and the solid line is a fit using the Wegner expansion 
and the rectilinear diameter law with the critical parameters given in Table HTTl 



parameters. 

Consider to the rectilinear diameter law 

l(p+ + p_] = Poo + bo(Too-T], (4.9) 

which is a empirical fit ^ extensively used for fluids in equilibrium. Here Poo 
and Tcx) denotes the critical density and the critical temperature, respectively. 
This, in principle, has no justification out of equilibrium. However, we found 
that our MC data nicely fit the diameters equation. We use this fact together 
with a universal scaling la-viH 

p+-p_ = ao(Too-T]P, (4.10) 

to accurately estimate the critical parameters. The simulation data in Fig. 14.111 
then yields the values in Table |4?T1 which are confirmed by the familiar log-log 
plots. Compared to the equilibrium critical temperature reported by Smit and 
Frenkel [55], one has that To /Too ~ 1-46, i.e., the change is opposite to the one 
for the DLG [TT] . This confirms the observation above that the field acts in the 
nonequilibrium LJ system favoring disorder. 

The fact that the order-parameter critical exponent is relatively small may 
already be guessed by noticing the extremely flat coexistence curve in Fig. 14.111 
This is similar to the corresponding curve for the equilibrium two-dimensional 
L J fluids [89l I117| 1118] , and it is fully consistent with the equilibrium Onsager 



^The first term of a Wegner- type expansion |116) . 
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Table 4.1: Critical parameters 



Poo 


Too 




0.321(5) 1 


0.314(1) 1 


0.10(8) 



value, (3 = 1/8. We therefore believe that our model belongs to the Ising uni- 
versality class. In any case, although the error bar is large, one may discard 
with confidence both the DLG value |3 w 1 /3 as well as the mean field value 
(3 = 1/2 — both cases would produce a hump visible to the naked eye in a plot 
such as the one in Fig. 14.111 One may argue that this result is counterintuitive, 
as our model apparently has the short-range interactions and symmetries that 
are believed to characterize the DLG. Further understanding for this difference 
will perhaps come from the statistical field theory. 

4.4 Conclusions 

In summary, the present (non-equilibrium) Lennard- Jones system, in which 
particles are subject to a constant driving field, has two main general features. 
On the one hand, this case is more convenient for computational purposes, than 
others such as, for instance, standard molecular-dynamics realizations of driven 
fluid systems. On the other hand, it seems to contain the necessary essen- 
tial (microscopic) physics to be useful as a prototypical model for anisotropic 
behavior in nature. 

This model is the natural extension to nonequilibrium anisotropic cases of 
the Lennard- Jones model which has played an important role in analysing equi- 
librium fluids. For zero field, it reduces to the familiar LJ fluid. Otherwise, 
it exhibits some arresting behavior, including currents and striped patterns. 
We have identified a process which seem to dominate early nucleation before 
anisotropic spinodal decomposition sets in. Interesting enough, they seem to be 
identical to the ones characterizing a similar situation in equilibrium. Surpris- 
ingly, we have also found that the model critical behavior is consistent with the 
Ising one for d = 2, i.e., (3 « 1 /8, but not with the critical behavior of the driven 
lattice gas, i.e., |3 ~ 1/3. This is puzzling. Moreover, this result is quantita- 
tively similar to the one in a related driven lattice system. Szolnoki and Szabo 
[91] showed that assuming nearest-neighbor repulsion under a driving field [93) 
and extending the dynamics to next-nearest-neighbor hops, the order-disorder 
transition which the model undergoes becomes second order for all drive, and 
belongs to the Ising universality class. Furthermore, it is straightforward to see 
that if one allows particles' spatial coordinates to vary continuously then the 
model reduces to a hard-disks gas under a driving field, and the second order 
phase transition vanishes. The transition becomes similar to the melting tran- 
sition for systems in equilibrium, which is described by the KTHNY scenario 
|119) . The fact that in both nonequilibrium diffusive models — which share sim- 
ilar dynamics — the critical behavior is consistent with Ising is rather surprising. 
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For instance, using bare arguments of statistical field theory, symmetries seem 
to bring our system closer to the nonequilibrium lattice model (the DLG) than 
to the corresponding equilibrium case (the LJ fluid). 

Therefore, a principal conclusion is that — regarding the modeling of complex 
systems — spatial discretization may change significantly not only morphologi- 
cal properties, but also critical properties. This is in contrast with the concept 
of universality in equilibrium systems, where critical properties are indepen- 
dent of dynamic details. The main reason for this disagreement might be the 
particle-hole symmetry violation in the driven Lennard-Jones fluid. However, 
to determine exactly this statement will require further study. 

In any case, the additional freedom of the present, off-lattice system, which 
in particular implies that the particle-hole symmetry which occurs in lattice 
models is violated — which induces the coexistence-curve asymmetry in Fig. 14.111 
in accordance with actual systems — are likely to matter more than suggested by 
some naive intuition. Indeed, the question of what are the most relevant ingre- 
dients and symmetries which determine unambiguously the universal properties 
in driven diffusive fluids is still open. Nevertheless, the above important dif- 
ference between the lattice and the off-lattice cases results most interesting as 
an unquestionable nonequilibrium effect; as it is well known, such microscopic 
detail is irrelevant to universality concerning equilibrium critical phenomena. 

Future research should address the development of a field-theoretical descrip- 
tion aiming at shedding light on the critical behavior of these models. Off-lattice 
models pose a new challenge for any available Langevin-type approach in which 
to test their viability. In principle, the anisotropic diffusive system approach 
described in the previous chapter appears to be more appropriate because it 
allows the consideration of the microscopic dynamical details. In fact, it is ex- 
pected that a careful analysis in the derivation process from the master equation 
to the final Langevin equation should rcvail the mechanisms that lead to Ising 
behavior. 

Further study of the present nonequilibrium L J system and its possible vari- 
ations is suggested. A principal issue to be investigated is the apparent fact 
that the full nonequilibrium situations of interest can be described by some 
rather straightforward extension of equilibrium theory. We here deal with some 
indications of this concerning early nucleation and properties of the coexistence 
curve. No doubt it would be interesting to compare more systematically the 
behavior of models against the varied phenomenology which was already re- 
ported for anisotropic fluids. This should also help a better understanding of 
nonequilibrium critical phenomena. 

Finally, one of our main aims at proposing this model has been motivate 
further research, both experimental and theoretical: we try to motivate new 
experiments on anisotropic decomposition besides theoretical work. In fact, 
our observations on the early-time separation process and structural properties 
are easily accessible by micro calorimetric and spectroscopy experiments, for 
instance. 
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Introduction 



The second part of this thesis focuses on rapid granular flows, also referred as 
granular gases [120\ \121\ 11221 11231 1124j . They are granular systems fluidized 
by the application of an external driving, e.g., by mechanical vibration |39[ 144] . 
shear [125; . electrostatically |126| . or others means [CT 1127) . In these systems, 
collisions between particles are considered nearly instantaneous compared with 
the mean free time. In fact, it is assumed (by definition) that the flow is gov- 
erned by binary collisions — as in molecular gases — , whereas three-body (or 
more) collisions are neglected. Instead of moving in many-particle blocks, each 
particle moves freely and independently of even its nearest neighbors. Fur- 
thermore, granular gases exhibit many states and instabilities, e.g., convection 
[48l l49l IT28] ■ that characterize molecular gases. In spite of these similarities with 
classical molecular fluids, there are significant differences between both which 
make granular gases strikingly different. The origin of differences lies, as noted 
in the general introduction (Chapter [T|), in the inelastic character of collisions 
between grains, and in the irrelevance of the thermal energy keT. 

One important coarse-grained (continuum) approach to modeling granular 
gases is hydrodynamics. Granular hydrodynamics has a great predictive power 
and looks ideally suitable for a description of large-scale — long-wavelength 
and long-time — patterns observed in rapid granular flows as, for instance, a 
plethora of clustering phenomena [43l EH [l29l [l30], vortices [3ll HnH [l32], os- 
cillons [39l I133j , and shocks |134j . This approach consists on a set of partial 
differential equations for the velocity v(r, t), temperature T(r, t) (average fluc- 
tuation kinetic energy), pressure p(r, t), and density (or number density n(r, t)) 
p(r, t) fields. Hydrodynamic equations are analogous to the Navier Stokes equa- 
tions for molecular fluids, but include a sink term in the energy equation which 
accounts for dissipation in collisions. For two-dimensional granular fluids, these 
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mass conservation, 
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V • p + ng momentum conservation, (5.1) 
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V • v^f + p : Vv - I 



energy conservation. 



Here, p is the stress tensor, w is the heat flux, I is the rate of energy loss by 
collisions per unit area per unit time, and g is a external force acting on the 
particles. Under some assumptions, more fundamental kinetic theory validates 
these equations and supplies constitutive relation|3 jl35| 11361 1137j . However, 
the exact criteria for the validity of (granular) kinetic theory and Navier-Stokes 
hydrodynamics [135[ 11381 1139] is still lacking. This is because first-principle 
derivations of universally applicable continuum theory of granular gas is not 
a simple task, even in the dilute limit. Unfortunately, the drastic simplifica- 
tion, provided by the binary collisions assumption, is far from sufficient for a 
derivation of hydrodynamics, and additional hypothesis are needed. Assuming 
a hard-disks system, and that the random motion of even neighboring particles 
are uncorrelated (molecular chaos assumption or stosszahlansatz) allows for 
derivations from the Boltzmann or Chapman-Enskog kinetic equation, properly 
generalized to account for the inelasticity of the particle collisions. An important 
additional assumption, made in the process of the derivation of hydrodynamics 
from the Chapman-Enskog equation, is scale separation |124) . Hydrodynamics 
demands (for inhomogeneous and/or unsteady) flows that the mean free path 
of the particles is much less than any characteristic length scale, and the mean 
free time between two consecutive collisions be much less than any characteris- 
tic time scale, aimed to be described hydrodynamically. All of these hypothesis 
are justified for moderate densities and for not too large inelasticity. These 
should be verified, in every specific system, after the hydrodynamic problem is 
solved and the characteristic length and time scales determined. In such case, 
systematic derivations of the constitutive relations — the equation of state and 
transport coefficients — of granular hydrodynamics from the Boltzmann or En- 
skog equation [1401 1141[ 1142] use an expansion in the Knudsen number and, in 
some versions of theory, in inelasticity [140) . 

A finite inelasticity immediately brings complications |124j . Correlations 
between particles, developing already at a moderate inelasticity, may invali- 
date the molecular chaos hypothesis |143| . The correlations become stronger 
as the inelasticity of the collisions increases. The normal stress difference and 
deviations of the particle velocity distribution from the Maxwell distribution 
also become important for inelastic collisions. As a result, the Navier-Stokes 
granular hydrodynamics should not be expected to be quantitatively accurate 
beyond the limit of small inelasticity. In any case, despite the severe limitations 
intrinsic to it, the nearly elastic case is conceptually important, just because 

^These are the explicit expressions of p, w, and I in terms of the basic fields n, T, and v. 
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hydrodynamics is supposed to give here a quantitatively accurate leading order 
theory. This is the regime (nearly elastic collisions) considered in the following 
Chapters |6] and [71 

A further set of difficulties for hydrodynamics arise at large densities. In 
mono-disperse systems strong correlations appear, already for elastic hard spheres, 
at the disorder-order transition. As multiple thermodynamic phases may coex- 
ist there |144) . a general continuum theory of hard sphere fluids, which has 
yet to be developed, must include, in addition to the hydrodynamic fields, an 
order-parameter field. Furthermore, even on a specified branch of the thermo- 
dynamic phase diagram, we do not have first-principle constitutive relations. 
Last but not least, kinetic theory of, for instance, a finite-density gas of elastic 
hard spheres in two dimensions is plagued by non-existence of transport coeffi- 
cients because of the long-time tail in the velocity pair autocorrelation function 
pl5l[T46] . 

Another potentially important limitation of the validity of granular hydro- 
dynamics (or, rather, of any continuum approach to rapid granular fiows) is due 
to the noise caused by the discrete nature of particles (discrete particle noise). 
This may play a dramatic role in rapid granular flow |1491 11501 1147] . 

Noise is stronger here than in classical molecular fluids simply because the 
number of particles is much smaller. In addition, noise can be amplified at 
thresholds of hydrodynamic instabilities as found, for example, in Rayleigh- 
Benard convection of classical fluids |148| or in granular gases |147| . A new 
challenge for theory is a quantitative account of this noise. A promising ap- 
proach at small or moderate densities is "Fluctuating Granular Hydrodynam- 
ics" : a Langevin-type theory that takes into account the discrete particle noise 
by adding a delta-correlated (gaussian) noise terms in the momentum and en- 
ergy equations, in the spirit of the Fluctuating Hydrodynamics by Landau and 
Lifshitz |15H I152j . At high densities, however, such a theory is again beyond 
reach. 

An alternative approach to granular flows, besides experiments, are simula- 
tions. Computer simulations have a valuable role in providing essentially exact 
results for certain problems which would otherwise be intractable. In addition, 
simulations help us to obtain the intuition we need in order to understand the 
complex behavior exhibited by granular material. Here (Chapters |6] and [71), we 
test the predictions from the hydrodynamic equations by comparing with Molec- 
ular Dynamics simulations [51] [551 1153] . Molecular Dynamics is a technique for 
computing the steady state and transport properties of a many-body system. 
This is in many respects very similar to real experiments. The general idea is 
simple: integrate numerically the microscopic equations of motion (Newton's 
equationll) simultaneously for all particles (1 < i < N) 

d^r- 

TUi =F|(ri,--- ,rN,vi,--- ,vn) (5.2) 



■^Typically, because of the macroscopic size of grains, quantum effects can be discarded 
with confidence. 
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This type of Molecular Dynamics is called force-based. However, there are sys- 
tems in which this scheme is very inefficient (although, of course, still applicable) 
|51j . In systems where most of the time each of particles propagates along a 
ballistic trajectory, the typical duration of a collision is much shorter than the 
mean time between successive collisions of a particle, i.e., binary collisions can 
be assumed. This is, in fact, what occurs in rapid granular flows. In such a 
case, if the post-collisional velocities are obtained from the pre-coUisional ve- 
locities a event-driven Molecular Dynamics simulation can be set uplf|. As a 
result, event-driven algorithms are much more efficient, in the sense of simu- 
lation speed, than force-based algorithms (integrating directly Eq. 15. 2p . We 
will employ event-driven molecular dynamic simulations to test the validity of 
hydrodynamic prediction in the following chapters. 

A standard prototypical model for driven granular gases, employed in many 
theories and simulations, as well as in experiments, is that of a monodisperse 
collection of uniform frictionless hard spheres (hard disks in two-dimensions) 
whose collisions are characterized by constant normal restitution coefficient M-jJj. 
Usually, this model is employed to study dense or diluted flows, with or without 
gravity, and with or without a interstitial fluid. Here the tangential degrees of 
freedom are neglected, i.e., =1, and only the normal restitution coefficient 

appears. However, it is worth mentioning that these two assumptions are 
ill-suited (see a nice discussion in [5T]). This simplification is owed exclusively 
to the huge theoretical complication which comes along with the incorporation 
of the correct (velocity dependent) tangential restitution coefficient. More com- 
plex models may involve nonspherical particles, in general, irregularly shaped 
particles, and are also polydisperse. In addition, except for astrophysical stud- 
ies [122) or experiments performed in vacuum [3 9) the grains are embedded in 
an ambient fluid. Nonetheless, the effect of a interstitial fluid, as air or wa- 
ter, can be sometimes neglected in determining many properties of the system 
|120) . This is particularly true when grain-grain interactions are stronger than 
fluid or grain-fluid interactions. But as the level of detail increases so does the 
computational effort. In spite of its shortcomings as a practical tool, the hard- 
spheres prototype enables the elucidation of many essential features which are 
the responsible for the observed macroscopic behavior. 

Regarding the mechanism which drives the particles, in Chapters |6] and [7] we 
consider rapidly vibrating boundaries. We will assume that the characteristic 
frequency of vibration of the wall is much larger than the collision rate of the 
grains. This guarantees the absence of correlations between successive collisions 
of particles with the vibrating wall. In addition, it will be assumed that the 
amplitude of the vibration is much smaller than the mean free path of the par- 
ticles, so collective motions in the system are avoided. Therefore, the vibrated 

There are previous works in which event-driven algorithms has been even appUed to dense 
granular systems |130) . 

■^The restitution coefficient typically depends weakly on the velocity, but using a constant 
simplifies the theoretical derivations greatly I51| . 



57 



wall can be represented as a hydrodynamic heat flux, i.e., as a "thermar wall. 
From a computational point of view, the thermal wall is implemented as follows: 
whenever a particle collides with the wall, it forgets its velocity and chooses a 
new one from a Maxwell distribution according to a certain, previously pre- 
scribed temperature T. It is important to notice that in order to avoid that the 
gas approaches the incorrect final temperature, the random normal velocities 
have to be chosen according to the probability distribution 



2T 



and the tangential velocities have to be determined from 
N(vt) - 



m 



exp 



rav^ 



(5.3) 



(5.4) 



See e.g. Ref. [51], pages 173-177, for detail. 



The first simulation of a system of particles colliding ineslatically and driven 
by a side themal wall was considered, in one dimension, by Du, Li, and Kadanoff 
|154| . For typical initial conditions, the system in Ref. [154 evolves to a state 
where the particles are separated into two groups. Almost all particles form a 
low-temperature cluster in a small region, while a few remaining particles move 
with high velocities. Clearly, this steady state cannot be described by granular 
hydrodynamics. The results in Ref. jl54| brought into question the validity of 
granular hydrodynamics in general. Nevertheless, some results [1551 1156j have 
showed that this "anomaly" seems to vanish in higher dimensions (our studies 
developed in Chapters [5] and [7] concern two-dimensional systems). Therefore, 
the validity of hydrodynamics merits further studies. 
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Chapter 6 



Phase separation of a 
driven granular gas in 
annular geometry 



In this chapter wc address phase separation of a monodisperse gas of inelas- 
tically coUiding hard disks confined in a two-dimensional annulus, the inner 
circle of which represents a "thermal wall" . When described by granular hy- 
drodynamic equations, the basic steady state of this system is an azimuthally 
symmetric state of increased particle density at the exterior circle of the annu- 
lus. When the inelastic energy loss is sufficiently large, hydrodynamics predicts 
spontaneous symmetry breaking of the annular state, analogous to the van der 
Waals-like phase separation phenomenon previously found in a driven granular 
gas in rectangular geometry. At a fixed aspect ratio of the annulus, the phase 
separation involves a "spinodal interval" of particle area fractions, where the gas 
has negative compressibility in the azimuthal direction. The heat conduction in 
the azimuthal direction tends to suppress the instability, as corroborated by a 
marginal stability analysis of the basic steady state with respect to small per- 
turbations. To test and complement our theoretical predictions we performed 
event-driven molecular dynamics simulations of this system. We clearly identify 
the transition to phase separated states in the simulations, despite large fluctua- 
tions present, by measuring the probability distribution of the amplitude of the 
fundamental Fourier mode of the azimuthal spectrum of the particle density. 
We find that the instability region, predicted from hydrodynamics, is always 
located within the phase separation region observed in the simulations. This 
implies the presence of a binodal (coexistence) region, where the annular state 
is metastable. 
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6.1 Introduction 

This chapter deal with a simple model of rapid granular flows, also referred 
to as granular gases: large assemblies of inelastically colliding hard spheres 
[Uni HH [HI [m [m IISI]- in the simplest version of this model the only 
dissipative effect taken into account is a reduction in the relative normal velocity 
of the two colliding particles, modeled by the coefficient of normal restitution, 
see below. 

As we explained in the previous chapter, under some assumptions a hydro- 
dynamic description of granular gases becomes possible. In short, the molec- 
ular chaos assumption allows for a description in terms of the Boltzmann or 
Enskog equations, properly generalized to account for the inelasticity of par- 
ticle collisions, followed by a systematic derivation of hydrodynamic equations 
|140[ I141[ 1142] , and for inhomogeneous flows hydrodynamics demands also scale 
separation. The implications of these conditions can be usually seen only a 
posteriori, after the hydrodynamic problem in question is solved, and the hy- 
drodynamic length/time scales are determined. 

We will restrict ourselves in this chapter (and also in Chapter [7]) to nearly 
elastic collisions and moderate gas densities where, based on previous stud- 
ies, hydrodynamics is expected to be an accurate leading order theory. Though 
quite restrictive, these assumptions allow for a detailed quantitative study (and, 
quite often, prediction) of a variety of pattern formation phenomena in granular 
gases. One of these phenomena is the phase separation instability |128[ I147[ 
I158| I159( I160| I16H I162j that was first predicted from hydrodynamics and then 
observed in molecular dynamic simulations. This instability arises already in a 
very simple, indeed prototypical setting: a monodisperse granular gas at zero 
gravity confined in a rectangular box, one of the walls of which is a "thermal" 
wall. The basic state of this system is the stripe state. In the hydrodynamic 
language it represents a laterally uniform stripe of increased particle density at 
the wall opposite to the driving wall. The stripe state was observed in experi- 
ment [129j . and this and similar settings have served for testing the validity of 
quantitative modeling |154[ 11551 1156] . It turns out that (i) within a "spinodal" 
interval of area fractions and (ii) if the system is sufficiently wide in the lateral 
direction, the stripe state is unstable with respect to small density perturbations 
in the lateral direction |128[ I160| I161j . Within a broader "binodal" (or coexis- 
tence) interval the stripe state is stable to small perturbations, but unstable to 
sufficiently large ones |158( 1162) . In both cases the stripe gives way, usually via 
a coarsening process, to coexistence of dense and dilute regions of the granu- 
late (granular "droplets" and "bubbles" ) along the wall opposite to the driving 
wall |158[ 11591 1162j . This far-from-equilibrium phase separation phenomenon is 
strikingly similar to a gas-liquid transition as described by the classical van der 
Waals model, except for large fluctuations observed in a broad region of aspect 
ratios around the instability threshold |147) . The large fluctuations have not 
yet received a theoretical explanation. 

This chapter addresses a phase separation process in a different geometry. 
We will deal here with an assembly of hard disks at zero gravity, colliding inelas- 
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tically inside a two-dimensional annulus. The interior wall of the annulus drives 
the granulate into a non-equilibrium steady state with a (hydrodynamically) 
zero mean flow. Particle collisions with the exterior wall are assumed elastic. 
The basic steady state of this system, as predicted by hydrodynamics, is the 
annular state: an azimuthally symmetric state of increased particle density at 
the exterior wall. The phase separation instability manifests itself here in the 
appearance of dense clusters with broken azimuthal symmetry along the exterior 
wall. Our main objectives are to characterize the instability and compute the 
phase diagram by using granular hydrodynamics (or, more precisely, granular 
hydrostatics, see below) and event driven molecular dynamics simulations. By 
focusing on the annular geometry, we hope to motivate experimental studies 
of the granular phase separation which may be advantageous in this geometry. 
The annular setting avoids lateral side walls (with an unnecessary /unaccounted 
for energy loss of the particles). Furthermore, driving can be implemented here 
by a rapid rotation of the (slightly eccentric and possibly rough) interior circl^ll- 

We organized the chapter as follows. Section l6^ deals with a hydrodynamic 
description of the annular state of the gas, which is completely described by 
three dimensionless parameters: the grain area fraction, the inelastic heat loss 
parameter, and the aspect ratio. As we will be dealing only with states with a 
zero mean flow, we will call the respective equations hydrostatic. 

A marginal stability analysis predicts a spontaneous symmetry breaking of 
the annular state. We compute the marginal stability curves and compare them 
to the borders of the spinodal (negative compressibility) interval of the system. 
In Section 16.31 we report event-driven molecular dynamics simulations of this 
system and compare the simulation results with the hydrostatic theory. In 
Section [6.41 we summarize the main results. 

6.2 Particles in an annulus and granular hydro- 
statics 

6.2.1 The density equation 

Let N hard disks of diameter d and mass m = 1 move, at zero gravity, inside an 
annulus of aspect ratio CI — Roxt/Rint, where Rext is the exterior radius and Rint 
is the interior one. The disks undergo inelastic collisions with a constant coef- 
ficient of normal restitution For simplicity, we neglect the rotational degree 
of freedom of the particles. In each inter-particle collision the kinetic energy 
is continuously transferred into heat, while the momentum is conserved. The 
(driving) interior wall is modeled by a thermal wall kept at temperature To, 
whereas particle collisions with the exterior wall are considered elastic. The en- 
ergy transferred from the thermal wall to the granulate dissipates in the particle 
inelastic collisions, and we assume that the system reaches a (non-equilibrium) 

^It is worth notice that the posed setting resembles morphology and dynamical processes 
in planetary rings 11631 . where clustering, spontaneous symmetry breaking, and oscillatory 
instabilities, between many others, also may occur. 
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steady state with a zero mean flow. We restrict ourselves to the nearly elas- 
tic limit by assuming a restitution coefficient close to, but less than, unity: 
1 — |j. <C 1 . This allows us to safely use granular hydrodynamics |124j . For a 
zero mean flow steady state the continuity equation is obeyed trivially, while 
the momentum and energy equations yield two hydrostatic relations: 

V-w(r) + l = 0, p = const, (6.1) 

where w is the local heat flux, I is the energy loss term due to inelastic collisions, 
and P — P(n, T) is the gas pressure that depends on the number density Ti(r) 
and granular temperature T(r). Despite the ineslatic collisions, if the macro- 
scopic field gradients are not very large, it seems natural to postulate linear 
relation between fluxes and "thermodynamic" forces. Then, we adopt the clas- 
sical Fourier relation for the heat flux w(r) = — KVT(r) (where k is the thermal 
conductivity), omitting a density gradient term y VTi(r). In the dilute limit this 
term was derived in Ref. |141j . It can be neglected in the nearly elastic limit 
which is assumed throughout the second part of this thesis. 
The momentum and energy balance equations read 

V • [KVT(r)] = I, p = const, (6.2) 

To get a closed formulation, we need constitutive relations for p(n, T), K{n, T) 
and I(n, T). We will employ the widely used semi-empiric transport coefficients 
derived by Jenkins and Richman [1351 from the Enskog kinetic equation, gen- 
eralized to account for inelastic collision losses, for moderate densities: 
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(6.3) 



and the equation of state first proposed by Carnahan and Starling |164) 

p=nT(1+2G), (6.4) 

where G = v(1 — 7^)/(1 — v)^ and y — n (7Td^/4) is the solid fraction. Let us 
rescale the radial coordinate by Ri^t and introduce the rescaled inverse density 

Z(r, 9) — nc/n(r, 0), where — 2/ ^V3d^^ is the hexagonal close packing 
density. The rescaled radial coordinate r now changes between 1 and CI = 
Rcxt/Rint, the aspect ratio of the annulus. As in the previous work |128| . Eqs. 
(|6.2p . (|6.4p and (|6.3p can be transformed into a single equation for the inverse 
density Z(r): 

V • [TiZWZ] = AS(Z) , (6.5) 
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where 
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The dimensionless parameter A = (27t/3)(1 — |x) (Rint/d)^ is the hydrodynamic 



(6.7) 



inelastic loss parameter. The boundary conditions for Eq. (j6.5p are 

az(i,e)/ae = o and VnZ(a,e)=o, 



The first of these follows from the constancy of the temperature at the (thermal) 
interior wall which, in view of the constancy of the pressure in a steady state, 
becomes constancy of the density. The second condition demands a zero normal 
component of the heat flux at the elastic wall. Notice that neither the pressure 
nor the temperature at the thermal wall enter in Eqs. (|6.5P - (I6.7I) [16111128] . This 
is a consequence of the fact that the hard-core interactions between particles 
does not introduce a characteristic energy scale. The temperature only sets the 
temperature scale in the system and affects the pressure, which is constant in 
the steady state problem. Finally, working with a fixed number of particles, we 
demand the normalization condition 



In 



a 



Z"^ (r, 0)rdrde =7Tf (D^-r 



where 



N 



is the area fraction of the grains in the annulus. Equations (|6.5|) - (|6.8|) deter- 
mine all possible steady state density profiles, governed by three dimensionless 
parameters: f. A, and O. 



6.2.2 Annular state 



The simplest solution of the density equation ()6.5|) is azimuthally symmetric 
(6-independent): Z — z(r). Henceforth we refer to this basic state of the system 
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Figure 6.1: The normalized density profiles obtained from the molecular dy- 
namics simulations (the dots) and hydrostatics (the line) for O = 2, A = 81.09, 
and f — 0.356 (equivalently, za — 2.351 ). The simulations were carried out with 
N = 1250 particles, |_l — 0.92, and Rjnt = 22.0. Also shown is a typical snapshot 
of the system at the steady state as observed in the simulation. 



as the annular state. It is determined by the following equations: 



[rJ-lzjz']' = rAQ(z), z'iCl] = 0, and 
z-^dr = (Q2-1)f/2, 



(6.9) 



where the primes denote r-derivatives. In order to solve the second order equa- 
tion (j6.9p numerically, one can prescribe the inverse density at the elastic wall 
Za = Combined with the no-flux condition at r = 11, this condition de- 

fine a Cauchy problem for z(t] |1591 1128) . Solving the Cauchy problem, one can 
compute the respective value of f from the normalization condition in Eq. (j6.9p . 
At fixed A and CI, there is a one-to-one relation between Zq and f. There- 
fore, an alternative parameterization of the annular state is given by the scaled 
numbers zq, A, and CI. The same is true for the marginal stability analysis 
performed in the next subsection. 

Figure 16.11 depicts an example of annular state that we found numerically. 
One can see that the gas density increases with the radial coordinat^, as ex- 
pected from the temperature decrease via inelastic losses (the velocities of the 
particles decrease), combined with the constancy of the (hydrodynamic) pres- 
sure throughout the system. The hydrodynamic density profile agrees well with 
the one found in our simulations, see below. 



2 See Refs. [42] US |T59l HSl] for clustering instability. 
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Figure 6.2: The main graph: the marginal stability curves k = k(f] (where 
k is an integer) for O = 1.5 and A = IC* (circles), A = 5 x IC* (squares), 
and A = 10^ (triangles). For a given A the annular state is stable above the 
respective curve and unstable below it, as indicated for A = IC*. As A increases 
the marginal stability interval shrinks. The lines connecting dots serve to guide 
the eye. The inset: the dependence of kn^ax on A^/^. The straight line shows 
that, at large A, k^ax cx A^/^. 

6.2.3 Phase separation 

Mathematically, phase separation manifests itself in the existence of additional 
solutions to Eqs. (|6.5p - (|6.8p in some region of the parameter space f , A, and CI. 
These additional solutions are not azimuthally symmetric^. Solving Eqs. (|6.5I) - 
(|6.8p for fully two-dimensional solutions is not easy |161j . One class of such 
solutions, however, bifurcate continuously from the annular state, so they can 
be found by hnearizing Eq. (j6.5p . as in rectangular geometry |161[ 1128] . In 
the framework of a time-dependent hydrodynamic formulation, this analysis 
corresponds to a marginal stability analysis which involves a small perturbation 
to the annular state. For a single azimuthal mode ~ sin(k9) (where k is integer) 
we can write Z(t, 9) = z(t) -|- e E(r)sin(k9), where E[r) is a smooth function, 
and £ <C 1 a small parameter. Substituting this into Eq. (|6.5|) and linearizing 
the resulting equation yields a k-dependent second order differential equation 
for the function r(r) = J^[Z(r)] E(t): 




(6.10) 



This equation is complemented by the boundary conditions 



r(i]=o and r'(a) = o. 



(6.11) 



•^As we shall show in Section 16.31 such instability also is observed in simulations. 



66 Phase separation of a driven granular gas in annular geometry 




0.2 0.4 0.6 0.8 1 1.2 1.4 

f x(a"-i)/2 

Figure 6.3: Two-dimensional projections on the A-f plane of the phase diagram 
at n = 1.5 (the solid line), = 3 (the dotted line), and O. = 5 (the dashed 
line). The inset shows more clearly, for = 3, that the marginal stability curve 
(the solid line) lies within the negative compressibility region (depicted by the 
dashed line). 



For fixed values of the scaled parameters f. A, and CI, Eqs. (j6.10|) and (|6.1ip 
determine a linear eigenvalue problem for k. Solving this eigenvalue problem 
numerically, one obtains the marginal stability hypersurface k = k(f, A, CI). For 
fixed A and CI, we obtain a marginal stability curve k = k(f). Examples of 
such curves, for a fixed CI and three different A are shown in Fig. 16.21 Each 
k = k(f) curve has a maximum kmax, so that a density modulation with the 
azimuthal wavenumber larger than kmax is stable. As expected, the instability 
interval is the largest for the fundamental mode k = 1 . The inset in Fig. 16.21 
shows the dependence of k^ax on A^^^. The straight line shows that, at large 
A, kmax oa A^/^, as in rectangular geometry |128j . 

Two-dimensional projections of the (f. A, Il)-phase diagram at three dif- 
ferent CI are shown in Fig. 16.31 for the fundamental mode. The annular state 
is unstable in the region bounded by the marginal stability curve and stable 
elsewhere. Therefore, the marginal stability analysis predicts loss of stability 
of the annular state within a finite interval of f, that is at frain(A, O) < f < 

fmax{A,a). 

The physical mechanism of this phase separation instability is the negative 
compressibility of the granular gas in the azimuthal direction, caused by the 
inelastic energy loss. To clarify this point, let us compute the pressure of the 
annular state, given by Eq. (j6.4l) . First we introduce a rescaled pressure P = 
p/(ni;To) and, in view of the pressure constancy in the annular state, compute 
it at the thermal wall, where T = To is prescribed and z[^] is known from our 
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numerical solution for the annular state. We obtain 



p(f,A,a) 



1 +2g(z(1)) 
z(1) 



The spinodal (negative compressibility) region is determined by the necessary 
condition for the instability: (dP/df)^^ n ^' whereas the borders of the spin- 
odal region are defined by (9P/9f)^Q — 0. That is, computing the deriva- 
tive we arrive at aP/9f = (9P/9z(1 ))(9z(1 )/9z(a))(9z(a)/9f). One can easily 
check that the first and third multipliers on the righthand side of this relation 
are always negative, thus the sign of (9P/9f)^ ^ is determined by the sign of 



Typical P(f) curves for a fixed 11 and several different A are shown in Fig. l6.4l 
One can see that, at sufficiently large A, the rescaled pressure P goes down with 
an increase of f at an interval fi < f < fa- That is, the effective compressibility 
of the gas with respect to a redistribution of the material in the azimuthal 
direction is negative on this interval of area fractions. By joining the spinodal 
points f] and (separately) fi at different A, we can draw the spinodal line 
for a fixed CI. As A goes down, the spinodal interval shrinks and eventually 
becomes a point at a critical point (Pc,fc], or (Ac,fc) (where all the critical 
values are H-dependent). For A < A^ P(f) monotonically increases and there 
is no instability. 

What is the relation between the spinodal interval (f i , fi) and the marginal 
stability interval (fmini fmax)? These intervals would coincide were the az- 
imuthal wavelength of the perturbation infinite (or, equivalently, k — > 0), so 
that the azimuthal heat conduction would vanish. Of course, this is not possi- 
ble in annular geometry, where k > 1 . As a result, the negative compressibility 
interval must include in itself the marginal stability interval (fmin.) fmax)- This 
is what our calculations indeed show, see the inset of Fig. 16.31 That is, a neg- 
ative compressibility is necessary, but not sufficient, for instability, similarly to 
what was found in rectangular geometry [128j . 

Importantly, the instability region of the parameter space is by no means 
not the whole region the region where phase separation is expected to occur. 
Indeed, in analogy to what happens in rectangular geometry |158[ 1162] . phase 
separation is also expected in a hinodal (or coexistence) region of the area frac- 
tions, where the annular state is stable to small perturbations, but unstable to 
sufficiently large ones. The whole region of phase separation should be larger 
than the instability region, and it should of course include the instability re- 
gion. Though we did not attempt to determine the binodal region of the system 
from the hydrostatic equations (this task has not been accomplished yet even 
for rectangular geometry, except in the close vicinity of the critical point [I62j ). 
we determined the binodal region from our molecular dynamics simulations re- 
ported in the next section. 
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Figure 6.4: The scaled steady state granular pressure P versus the grain area 
fraction f for Q = 1.5 and A = 1.1 x 10^ (the dotted line), A = 1.5 x 10^ 
(the dash-dotted line), A = 3.5 x 10^ (dashed line), and A = 5 x 10'* (the solid 
line). The inset shows a zoom- in for A = 3.5 x 10^. The borders fi and fi 
of the spinodal interval are determined from the condition (9P/9f)^ n ~ ^■ 
The thick solid line encloses the spinodal balloon where the effective azimuthal 
compressibility of the gas is negative. 

6.3 Molecular Dynamics Simulations 
6.3.1 Method 

We performed a series of event-driven molecular dynamics simulations of this 
system using an algorithm described by Poschel and Scliwager [51 and sketched 
in the previous chapter. Simulations involved N hard disks of diameter d = 1 
and mass ra = 1 . After each collision of particle i with particle j , their relative 
velocity is updated according to 

v(. =vij-(1+^)(vtj-t|j)fij, (6.12) 

where Vij is the precoUisional relative velocity, and f ij = I'lj / |fij | is a unit vector 
connecting the centers of the two particles. Particle collisions with the exterior 
wall r = Rext are assumed elastic. The interior wall is kept at constant temper- 
ature To that we set to unity. This is implemented as explained in Chapter O 
When a particle collides with the wall it forgets its velocity and picks up a new 
one from a proper Maxwellian distribution with temperature To . The time scale 
is therefore d(ra/To)^/^ = 1. The initial condition is a uniform distribution 
of non-overlapping particles inside the annular box. Their initial velocities are 
taken randomly from a Maxwellian distribution at temperature To = 1 . In all 
simulations the coefficient of normal restitution \i — 0.92 and the interior radius 
Rint/d = 22.0 were fixed, whereas the the number of particles 527 < N < 7800 
and the aspect ratio 1 .5 < D. < 6 were varied. In terms of the three scaled 
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Figure 6.5: Typical steady state snapshots for N — 1250 and Q = 6 (a); N = 
5267 and a = 3 (b), and N = 6320 and O = 6 (c). Panels (a) and (b) 
correspond to annular states of the hydrostatic theory, whereas panel (c) shows 
a broken-symmetry (phase separated) state. 

hydrodynamic parameters the heat loss parameter A = 81 .09 was fixed whereas 
f and O varied. 

To compare the simulation results with predictions of our hydrostatic theory, 
all the measurements were performed once the system reached a steady state. 
This was monitored by the evolution of the total kinetic energy (1 /2) ^ 
which first decays and then, on the average, stays constant. 

6.3.2 Steady States 

Typical steady state snapshots of the system, observed in our molecular dynam- 
ics simulation, are displayed in Fig. 16.51 Panel (a) shows a dilute state where the 
radial density inhomogeneity, though actually present, is not visible by naked 
eye. Panels (b) and (c) do exhibit a pronounced radial density inhomogeneity. 
Apart from visible density fluctuations, panels (a) and (b) correspond to annu- 
lar states. Panel (c) depicts a broken-symmetry (phase separated) state. When 
an annular state is observed, its density profile agrees well with the solution of 
the hydrostatic equations ()6.5|) - (|6.8|) . A typical example of such a comparison 
is shown in Fig. 16.11 

Let us fix the aspect ratio CI of the annulus at not too a small value and 
vary the number of particles N. First, what happens on a qualitative level? The 
simulations show that, at small N, dilute annular states, similar to snapshot (a) 
in Fig. 16. 5[ are observed. As IM increases, broken-symmetric states start to 
appear. Well within the unstable region, found from hydrodynamics, a high 
density cluster appears, like the one shown in Fig. 16.5b . and performs an erratic 
motion along the exterior wall. As N is increased still further, well beyond the 
high-f branch of the unstable region, an annular state reappears, as in Fig. 16.5b . 
This time, however, the annular state is denser, while its local structure varies 
from a solid-like (with imperfections such as voids and line defects) to a liquid- 
like. 

To characterize the spatio-temporal behavior of the granulate at a steady 
state, we followed the position of the center of mass (COM) of the granulate. 
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Figure 6.6: Typical steady state snapshots (the left column; particle size looks 

larger than it is) and the temporal evolution of the COM (the middle column) 
and of the squared amplitude of the fundamental Fourier mode (the right col- 
umn). The temporal data are sampled each 150 collisions per particle. Each 
row corresponds to one simulation with the indicated parameters. The vertical 
scale of panels (a) and (b) was stretched for clarity. 
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Several examples of the COM trajectories are displayed in Fig. 16.61 Here cases 
(a) and (b) correspond, in the hydrodynamic language, to annular states. There 
are, however, significant fluctuations of the COM around the center of the an- 
nulus. These fluctuations are of course not accounted for by hydrodynamic 
theory. In case (b), where the dense cluster develops, the fluctuations are much 
weaker that in case (a). More interesting are cases (c) and (d). They corre- 
spond to broken-symmetry states: well within the phase separation region of 
the parameter space (case c) and close to the phase separation border (case 
d). The COM trajectory in case (c) shows that the granular "droplet" per- 
forms random motion in the azimuthal direction, staying close to the exterior 
wall. This is in contrast with case (d) , where fluctuations are strong both in the 
azimuthal and in the radial directions. Following the actual snapshots of the 
simulation, one observes here a very complicated motion of the "droplet", as 
well as its dissolution into more "droplets" , mergers of the droplets etc. There- 
fore, as in the case of granular phase separation in rectangular geometry |147| . 
the granular phase separation in annular geometry is accompanied by consider- 
able spatio-temporal fluctuations. In this situation a clear distinction between 
a phase-separated state and an annular state, and a comparison between the 
simulations and hydrodynamic theory, demand proper diagnostics. We found 
that such diagnostics are provided by the azimuthal spectrum of the particle 
density and its probability distribution. 



6.3.3 Azimuthal Density Spectrum 



Let us consider the (time-dependent) rescaled density field v(r, 9, t) = n(r, 9, t]/nc 
(where r is rescaled to the interior wall radius as before), and introduce the in- 
tegrated field V(9,t): 



^(9,t) = 



v(r, 9,t)rdr. 



In a system of N particles, V(9,t] is normalized so that 



■27t 



'^(9,t) d9 



N 



(6.13) 



(6.14) 



Because of the periodicity in 9 the function V(9, t) can be expanded in a Fourier 
series: 



V(9,t) = ao + ^ [Qk(t)cos(k9) +b,^(t)sin(k9)] 



(6.15) 



k=l 



where qq is independent of time because of the normalization condition (16.141) . 
We will work with the quantities 



Ai{t) = ai[t)+hi(t], k>1. 



(6.16) 



For the (deterministic) annular state one has Aj^ = for all k > 1 , while for a 
symmetry-broken state > 0. The relative quantities A^(t)/QQ can serve as 
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Figure 6.7: The normahzed probabihty distribution functions P] (A^/ag) for 
0=3 (the left column) and D = 5 (the right column) at a different number of 
particles. 
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Table 6.1: The averaged squared relative amplitudes (A^(t))/ao for the first 
three modes k = 1,2 and 3. (a) N = 2634, O = 3, (b) N = 5267, 0=4, (c) 
N = 1000, a = 2.25, and (d) N = 1250, n = 3. 



k 


(a) 


(b) 


(c) 


(d) 


1 


0.66 ± 0.05 


0.39 ± 0.04 


0.30 ± 0.08 


0.77 ± 0.05 


2 


0.04 ± 0.02 


0.05 ± 0.02 


0.07 ±0.01 


0.28 ± 0.09 


3 


0.03 ± 0.02 


0.03 ± 0.03 


0.02 ± 0.02 


0.11 ±0.08 



measures of the azimuthal symmetry breaking. As is shown in Table IFTTl A^(t) 
is usually much larger (on the average) that the rest of A^(t). Therefore, the 
quantity A^(t]/aQ is sufficient for our purposes. 

Once the system relaxed to a steady state, we followed the temporal evolution 
of the quantity A^ /qq . Typical results are shown in the right column of Fig. 16.61 
One observes that, for annular states, this quantity is usually small, as is the 
cases (a) and (b) in Fig. 16.61 For broken-symmetry states A^ is larger, and it 
increases as one moves deeper into the phase separation region. (Notice that 
the averaged value of A^/qq in (c) is larger than in (d), which means that (c) is 
deeper in the phase separation region.) Another characteristics of A^(t)/aQ is 
the magnitude of fluctuations. One can notice that, in the vicinity of the phase 
separation border the fluctuations are stronger (as in case (d) in Fig. 16. 6[) . 

All these properties are encoded in the probability distribution Pi of the 
values of (Ai/ao)^: the ultimate tool of our diagnostics. Figure shows 
two series of measurements of this quantity at different N: for CI — 3 and 
CI — 5. By following the position of the maximum of Pi we were able to to 
sharply discriminate between the annular states and phase separated states and 
therefore to locate the phase separation border. When the maximum of Pi 
occurs at the zero value of (Ai/ao)^ (as in cases (a) and (d) and, respectively, 
(e) and (h) in Fig. 16.71) . an annular state is observed. On the contrary, when 
the maximum of Pi occurs at a non-zero value of (Ai/qq) (as in cases (b) 
and (c) and, respectively, (f) and (g) in Fig. 16. 7|) . a phase separated state is 
observed. In each case, the width of the probability distribution (measured, 
for example, at the half-maximum) yields a direct measure of the magnitude 
of fluctuations. Near the phase separation border, strong fluctuations (that is, 
broader distributions) are observed, as in case (c) of Fig. 16.71 

Using the position of the maximum of Pi as a criterion for phase separation, 
we show, in Fig. 16. 8[ the Cl — i diagram obtained from the molecular dynamics 
simulations. The same figure also depicts the hydrostatic prediction of the 
instability region. One can see that the instability region is located within the 
phase separation region, as expected. 
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Figure 6.8: The fl-f phase diagram for A = 81.09. The solid curve is given by 
the granular hydrostatics: it shows the borders of the region where the annular 
state is unstable with respect to small perturbations. The filled symbols depict 
the parameters in which phase separated states are observed, whereas the hollow 
symbols show the parameters at which annular states are observed. The dashed 
line is an estimated binodal line of the system. 

6.4 Conclusions 

We combined equations of granular hydrostatics and event-driven molecular dy- 
namics simulations to investigate spontaneous phase separation of a monodis- 
perse gas of inelastically colliding hard disks in a two-dimensional annulus, the 
inner circle of which serves as a "thermal wall" . A marginal stability analysis 
yields a region of the parameter space where the annular state — the basic, az- 
imuthally symmetric steady state of the system — is unstable with respect to 
small perturbations which break the azimuthal symmetry. The physics of the 
instability is negative effective compressibility of the gas in the azimuthal di- 
rection, which results from the inelastic energy loss. Simulations of this system 
show phase separation, but it is masked by large spatio-temporal fluctuations. 
By measuring the probability distribution of the amplitude of the fundamen- 
tal Fourier mode of the azimuthal spectrum of the particle density we were 
able to clearly identify the transition to phase separated states in the simula- 
tions. We found that the instability region of the parameter space, predicted 
from hydrostatics, is located within the phase separation region observed in the 
molecular dynamics simulations. This implies the presence of a binodal (coex- 
istence) region, where the annular state is metastahle, similar to what has been 
found in rectangular geometry [158| 1162] . We hope our results will stimulate 
experimental work on the phase separation instability. 

Finally, we have also investigated an alternative setting in which the exterior 
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Figure 6.9: Tlie marginal stability lines for our main setting (the solid line) and 
for an alternative setting in which the thermal wall is at r = Rgxt and the elastic 
wall is at r = Ri„f (the dashed line). 

wall is the driving wall, while the interior wall is elastic. The corresponding 
hydrostatic problem is determined by the same three scaled parameters f , A and 
fl, but the boundary conditions must be changed accordingly. Here azimuthally 
symmetric clusters appear near the (elastic) interior wall. Symmetry breaking 
instability occurs here as well. We found very similar marginal stability curves 
here, but they are narrower (as shown in Fig. 16. 9p than those obtained for our 
main setting. 



76 Phase separation of a driven granulai' gas in annular geometry 



Chapter 7 

Close— packed granular 
clusters 

Dense granular clusters often behave like macro-particles. We address this inter- 
esting phenomenon in a model system of inelastically colliding hard disks inside 
a circular box, driven by a thermal wall at zero gravity. Molecular dynamics 
simulations show a close-packed cluster of an almost circular shape, weakly fluc- 
tuating in space and isolated from the driving wall by a low-density gas. The 
density profile of the system agrees well with the azimuthally symmetric solution 
of granular hydrostatic equations employing constitutive relations by Grossman 
et al, whereas the widely used Enskog-type constitutive relations show poor 
accuracy. We find that fluctuations of the center of mass of the system are 
Gaussian. This suggests an effective Langevin description in terms of a macro- 
particle, confined by a harmonic potential and driven by a delta-correlated noise. 
Surprisingly, the fiuctuations persist when increasing the number of particles in 
the system. 

7.1 Introduction 

This chapter addresses granular hydrodynamics and fiuctuations in a simple 
two-dimensional granular system under conditions when existing hydrodynamic 
descriptions [1401 1141[ 1142] break down because of large density, not large in- 
elasticity. In view of the difficulties mentioned in Chapters [5] and |6l attempts of 
a first-principle description of hydrodynamics and fluctuations should give way 
here to more practical, empiric or semi-empiric, approaches. One such approach 
to a hydrodynamic (or, rather, hydrostatic) description was suggested by Gross- 
man et al. |156) . In the present chapter, we put it into a test in an extreme 
case when macro-particles (granular clusters with the maximum density close 
to the hexagonal close packing) form. The model system we are dealing with 
was first introduced by Esipov and Poschel |155j . It is an assembly of N ^ 1 
identical disks of mass m, diameter d and coefficient of normal restitution \i, 
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placed inside a circular box of radius R at zero gravity (typical steady state 
configurations are shown in Fig. I7.1|) . The circular wall of the box, which sup- 
ply energy continuously to the granulate, is kept at constant temperature To, 
i.e., the vibrating boundary is represented as a thermal wall (under conditions 
specified in Chapter [5]). There are not additional external forces acting on the 
system. 

We used a model of inelastically colliding hard disks to develop kinetic and 
hydrodynamic descriptions. We measure, using event-driven molecular dynam- 
ics simulations, the radial density profiles of the system, including the close- 
packed part. Furthermore, we solve numerically a set of granular hydrostatic 
equations which employ the constitutive relations by Grossman et al. and show 
that, in a wide range of parameters, there is good agreement between the two. 
A marginal stability analysis will show that there are no steady-state solutions 
with broken azimuthal symmetry which would bifurcate from the azimuthally 
symmetric solution. We also show that, for the same setting, the Enskog-type 
CRs [135] perform poorly. Finally, we investigate in some detail, for the first 
time, fluctuations of the macro-particles, by measuring the radial probability 
distribution function (PDF) of the center of mass of the system. These fluctua- 
tions turn out to be Gaussian, which suggests an effective Langevin description 
of the system in terms of a macro-particle, confined by a harmonic potential 
and driven by a white noise. Surprisingly, the fiuctuations persist as the number 
of particles in the system is increased. 



7.2 Event— driven molecular dynamics simulations 

As we outlined in Chapter [21 we employed a standard event-driven molecu- 
lar dynamics algorithm ^51]. One of the simplest models one can employ to 
study fiows dominated by binary particle collisions |124[ I157j consists of N in- 
elastically colliding hard disks of mass ra and diameter d. Rotational degree 
of freedrom is neglected, and thus inelasticity is modeled by the coefficient of 
normal restitution |j. < 1, which is considered constant (velocity independent). 
When two particles collide (say 1 < i, j < N), their tangential velocities are un- 
changed, whereas the relative normal velocity is decreased. Using momentum 
conservation, the post-collisional velocities in terms of the pre-coUisional read 

(7-1) 

where primed quantities stand for post-collisional velocities, and 
tij = (fi — rj / |ri — f j I is unitary vector oriented along the direction connecting 
their centerf[j. The circular boundary supply continuously energy to the system 
as if were driven by a vibrating wall. If the frequency of the oscillating wall 



Notice that Eq. I|7.ip is equivalent to Eq. Ij6.12|l . but expressed in a more detailed way. 
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Figure 7.1: Snapshots of the system (top, from left to right) for N — 1716, 
2761 and N — 5055. Also shown are magnifications (bottom) of the indicated 
areas in the rightmost snapshot. Notice the coexistence of the different kinds of 
packings. 

is large enough, the energy supply may be modelled by a heated wall. The 
thermal wall implementation |51) is simple: whenever a particle touches the 
wall its velocity is chosen from a Maxwell distribution according to To (which 
is measure in units of energy), that we set to unity. We put m — 1 and fixed 
R/d — 100 and \i — 0.888, while the total number of particles N served as the 
control parameter in the simulations. So the scaled time unit is d (ra/To) ^ ^'^ — ^ . 

We were mostly interested in a hydrodynamic (low Knudsen number) regime, 
when the mean free path is small compared to the system size and the mean free 
time is small compared with any hydrodynamic time scale. This requires N 3> 
R/d = 100. For N in the range of a few hundred, we observed a dilute granular 
gas with an increased density in the center of the box (although not observable to 
naked eye). The clustering in the center becomes more pronounced as N grows. 
The clustering can be easily explained in the hydrodynamics language: because 
of the inelastic collisions the granular temperature goes down as one moves away 
from the wall toward the center of the box. This, combined with the constancy 
of the pressure throughout the system, causes an increased particle density at 
the center. As N increases further, the particle density in the center approaches 
the hexagonal close packing value tie = 2/(-\/3d^). Figure [TT] shows snapshots 
of the system for three different, but sufficiently large, values of N. A perfect 
hexagonal packing is apparent nearby the center. Movies of these simulations 
show that the macro-particle (close-packed cluster) position fluctuates around 
the center of the box, while the cluster shape fluctuates around a circular shape. 
Notice that both kind of fluctuations, in particular the off-center location if the 
cluster, can not be explained hydrodynamically. 
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Our diagnostics was focused on two radial distributions: the number density 
of the particles n(r) and the PDF of a radial position of the center of mass 
of the system P(rTa), defined below. To measure n(r) (likewise in Chapter |6]), 
we recorded all particle coordinates in intervals of 100,000 particle collisions. 
For each snapshot we introduced bins in the form of concentric circular rings 
of width 0.5, centered in the center of mass of the system. For each particle 
we determined the fraction of its area falling in each of the rings, summed the 
contributions of all particles to each bin, and divided the result by the bin area. 
The resulting radial profile was averaged over many snapshots (typically 1000). 
Similar method is employed to measure T'[rm). As we are interested in steady- 
state distributions, we disregarded initial transients. This was monitored by the 
time-dependence of the average kinetic energy of the particles (the first decayed 
and then approached an almost constant value) and the time-dependence of the 
center of mass itself. 



7.3 Hydrostatic theory 

As the fluctuations are relatively weak, it is natural to start with a purely 
hydrodynamic description. For a zero mean flow this is a hydrostatic theory, 
as it operates only with (time-independent) granular density ■n.(r), temperature 
T(r) and pressure p(r). The energy input at the thermal wall is balanced by 
dissipation due to inter-particle collisions, and one can employ the momentum 
and energy balance equations: 

p — const , 

^ 7.2 
V-{kVT)=I. ^ ' 

Here k is the thermal conductivity and I is the rate of energy loss by collisions. 
Note that the heat flux, entering the thermal balance in Eq. (|7.2p . does not 
include an inelastic term, proportional to the density gradient [1401 11411 1142j . 
In the nearly elastic limit 1 — |j. <C 1 , that we are interested in, this term can 
be neglected. The boundary condition at the thermal wall is T(r = R, 9) = To, 
where r and 6 are polar coordinates with the origin at the center of the box. 

To proceed from here and compute numerical factors, we need constitutive 
relations: an equation of state p = p (n., T) and relations for k and I in terms 
of n and T. As we attempt to describe close-packed clusters, the first-principle 
standard techniques, based on the Boltzmann or Enskog equations, are inap- 
plicable. Grossman et al. [156] suggested a set of semi-empiric relations in 
two dimensions, that are valid for all densities, all the way to hexagonal close 
packing. Their approach ignores possible coexistence beyond the disorder-order 
transition and assumes that the whole system is on the thermodynamic branch 
extending to the hexagonal close packing. The delicate issue of phase coexis- 
tence (liquid-like phase, close-packed phase with multiple domains) occur here 
is close analogy to the system of elastic hard spheres (1221 11441 1166] . Gross- 
man et al. employed free volume arguments in the vicinity of the close packing, 
and suggested an interpolation between the hexagonal-packing limit and the 
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well-known low-density relations. The resulting global equation of state and 
constitutive relations [156j read 



rir 



n 



p = nT 

n-c — n 

(jn(al + d)2Ti/2 



(7.3) 
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Here I is the mean free path, which is given by an interpolation formula [156] 
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(7.4) 



and a = 1 — (3/8)^/^. These relations include three dimensionless numerical 
factors of order unity: a, y and cr, where the latter drops out from the steady- 
state problem. They can be calculated exactly from the velocity distribution 
function. However, this distribution is poorly understood |124) . Grossman et 
al. determined the optimum values a = 1.15 and y — 2.26, verified by a 
detailed comparison with molecular dynamics simulations of a system of inelastic 
hard disks in a rectangular box without gravity, driven by a thermal wall, and 
numerical solutions of the hydrostatic equations (j7.2[) in rectangular geometry. 
We adopted the same values of and y in our calculations for the circular 
geometry. 

Employing Eqs. (|7.3|) and (|7.4[) we can reduce Eqs. (|7.2[) to a single equation 
for the rescaled inverse density z(r, 6) = nc/n(r, 6). In the rescaled coordinates 
r/R — > r the circle's radius is 1 , and the governing equation becomes 
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(7.5) 



where 
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(7.6) 



and A = (32/3y) (R/d)^ (1 — ^l^) is the hydrodynamic inelasticity parameter, 
first introduced in Ref. |161| in the context of rectangular geometry As the 
total number of particles N is fixed, z^^ (r, 0) satisfies a normalization condition: 



In 



d0 



dr r z ^ (r, 0) = Ttf , 



(7.7) 



^The value of a differs in systems with gravity (a ~ 0.6) |130| . 

''Notice that, no matter how small (but finite) the inelasticity is, the parameter A goes to 
infinity in the thermodynamic limit R — > 00. 
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Figure 7.2: Scaled density n/ric versus scaled radius r/R as observed in simu- 
lations (circles) and predicted by the hydrostatic theory with the constitutive 
relations by Grossman et al. |156) (solid lines) and with the Enskog-type rela- 
tions [135) (dashed lines) for four different values of the total number of particles 
N. For the rest of parameters please see the text. The dotted lines indicate 
n/ric = 1 . 



7.3 Hydrostatic theory 



83 



where 



In VR 



(7.8) 



is the average area fraction of the particles. 

Most interesting state among the one-dimensional states is the azimuthally 
symmetric state (9-independent), which is indeed one of the possible solutions 
of Eq. (|7.2p . In this case, assuming azimuthal symmetry, we rewrite Eq. ()7.2[) 
as 



r dr 



TK n,T) — 
dr 



p(n, T) — const , 
= I(n,T). 



(7.9) 



Equations (|7.9p can be reduced to a single equation for the scaled inverse density 
z(t) = rLc/n(r). In the rescaled coordinate r, Eqs. (|7.9p reduces to 



r dr 



rJ^ z — 
dr 



(7.10) 



The fixed total number of particles yields a normalization condition, which 
read 



z-i(r)rdr = f/2. 



(7.11) 



Equations (|7.10|) and (j7.11|) . together with the obvious boundary condition 
dz/dT|r=o — (this is a consequence of the fact that the temperature gradi- 
ent vanish at the center of the system), form a complete set. The hydrostatic 
problem is completely determined by two scaled parameters A and f . This can 
be solved numerically as explained in Chapter |6] for the annular geometry. That 
is, in order to solve Eq. ()7.10p . instead of consider the area fraction f one can 
prescribe the inverse density at the center of the system zq . This condition, com- 
bined with the no-flux condition at r = define a Cauchy problem for z^^ (r) 
|128j . Solving the Cauchy problem one can compute the respective value for f 
from Eq. (|7.1ip . At fixed A, zo is a strictly monotonically decreasing function 
of f . Therefore, an alternative parameterization of the azimuthally symmetric 
state is given by the scaled numbers A and zq. The same considerations keep 
for the marginal stability analysis (see below) . 



In our event-driven simulations we varied the number of particles N . Accord- 
ingly, the scaled parameter f is varied, while A — 9980 was constant. Figure [72] 
shows a comparison of the hydrostatic radial density profiles with the radial 
density profiles obtained in simulations, for four different values of N. The sim- 
ulations profiles were averaged over 1000 uncorrelated configurations. It can be 
seen from Fig. 17.21 that for N — M^G the theory overestimates the density in 
the center of the box. This is expected as, for relatively small N, the maximum 
density is considerably less than the close-packing density, and the accuracy 
of the constitutive relations by Grossman et al. is not as good. For larger N 
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Figure 7.3: Maximum scaled density nmax/iT-c as a function of N as predicted 
by the hydrostatic theory and observed in simulations. 



the agreement rapidly improves. As Fig. 17.31 shows, the agreement between the 
hydrostatic theory and the simulations persists for the maximum radial den- 
sities at different number of particles N. The agreement is excellent for the 
densities approaching the close packing density, whereas for small N, there is 
small deviations. In absolute scale, however, the agreement is very good. We 
also computed the density profiles using another set of semi-empiric constitutive 
relationqj, those obtained in the spirit of Enskog theory |135| . These read. 
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(7.12) 



and the equation of state first proposed by Carnahan and Starling |164) 

p=nT(1+2G), (7.13) 
where G = v(l — -y^j/ll — v)^ and v = n (7td^/4) is the solid fraction. One 



can clearly see from Fig. 17.21 that the Enskog- type constitutive relations predict 
unphysically high densities in the cluster. 



Are the azimuthally symmetric states stable with respect to small perturba- 
tions? We performed marginal stability analysis to find out whether there are 
steady-state solutions with broken azimuthal symmetry, 2(r, 6) that bifurcate 
continuously from an azimuthally symmetric solution z(r). Under additional 
assumption that the possible instability of the azimuthally symmetric state is 
purely growing (that is, not oscillatory), the marginal stability analysis yields 



*These were already employed in Chapter |6] 
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the instabihty borders. The marginal stabihty analysis goes along the same 
lines as that developed for the rectangular geometry jl28( I162[ I16H I159( I147j . In 
the framework of time-dependent hydrodynamics, this corresponds to marginal 
stability analysis of the azimuthally symmetric solution with respect to small 
perturbations along the azimuthal coordinate 6. Let us search a steady-state 
close to an azymuthally-symmetric state: 

2(r, 9) = z{r) + e 4)k(r) sin(ke) , (7.14) 

As 2(r, 4) + 2n] — £(t, (})], k must be an integer which can be chosen to be 
non-negative. Substituting Eq. (I7.14p into Eq. (|7.5p and linearizing around the 
azimuthally symmetric state z(r) with respect to the small correction e ^ 1 , 
we obtain a linear eigenvalue problem, where k = k(f , A) plays the role of the 
eigenvalue: 

where Ck(i") — 4)k(r)J^(z]. This equation is complemented by the boundary 
conditions 

Ck(0)=0 Ck(1)=0, (7.16) 

and can be solved numerically. Let us ignore for a moment the quantization of 
the eigenvalue k and, while looking for k, assume that it is a (positive) real num- 
ber. In that case, a numerical solution yields, for a fixed A, a curve k — k(f), 
see Fig. 17.41 for two examples. At a fixed k, the azimuthally symmetric state is 
unstable within an interval of area fractions. The foots of one such curve corre- 
sponds to the (hypothetical) case when k tends to zero (that is, the azimuthal 
wavelength tends to infinity). The instability interval becomes narrower when 
k is increased, and it shrinks to a point at a maximum kraax, signalling that 
a density modulation with a sufficiently short azimuthal wavelength should be 
stable for all f. For example, when A = 5 x IC*, the marginal stability curve 
k — k(f ) has its maximum at kraax ~ 0-38 (see Fig. 17. 4p that is less than unity. 
Going back to the physical case, where k is a positive integer, we see that all the 
values for k, determined from the eigenvalue problem, are unphysical. That is, 
there are no solutions to the eigenvalue problem that would satisiy the bound- 
ary conditions and the quantization condition k = 1 , 2, 3 ... . We observed a 
similar behavior for values of A up to 10^. Though km ax increases with A, the 
increase is extremely slow: slower than logarithmical (see the inset of Fig. 17. 4|) . 
These numerical results strongly indicate that the azimuthally symmetric states 
are stable with respect to small perturbations. This is in marked contrast with 
the presence of bifurcating states with broken symmetry in similar settings of a 
granular gas driven by a thermal wall, but in rectangular [1281 11621 ll6H[159lll47j 
and annular (Chapter |6]) geometries. The absence of the bifurcating states with 
broken symmetry in the circular geometry gives a natural explanation to the 
persistence of circle-shaped cluster shapes as observed in our MD simulations. 
This prediction is very robust, as it is independent of the values a and y, and 
even of the constitutive relations employed. For instance, the same result is 
found by employing Enskog-type relations |135j (see Eq. (|7.12|) '). 
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Figure 7.4: The marginal stability curves k = k[f) for A = 5 x lO'* and A = 10^, 
as indicated. The shaded area denotes in both cases the (unphysical) instability 

region that is obtained if one ignores the quantization of k: k = 1 , 2, The 

inset shows kraax as a function of ln(A). 
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7.4 Macro— particle fluctuations 

Now we turn to stationary fluctuations, that is, after transients die out. To 
better characterize the fluctuation-dominated behavior of this system, we com- 
pute the radial coordinate rm.[t) of the center of mass of the system. The 
radial probability distribution function P(TTn.,t) is normalized by the condition 
271 Jq P(Tm^, t) r^dtiT^ = 1, where we have returned to the dimensional coordi- 
nate. Typical molecular dynamics results are presented in Fig. 17.51 which shows 
logP(rm^] versus (tm/R)^, for four different values of N. The observed straight 
lines clearly indicate a Gaussian distribution, both for small and large N values. 

This finding strongly suggests a Langevin description of the macro-particle. 
One can consider the macro-particle performing an over-damped motion in a 
(central-symmetric) confining harmonic potentialf] Ufrm.) = and driven 

by a (delta-correlated) discrete-particle noise r\ (t) . The Langevin equation for 
this problem reads 

Hf„,-f-kr^ =ri(t), (7.17) 

where h is the damping rate, (ri(t)ri(t')) — 2hr5(t — t'), and V is an effective 
magnitude of the discrete-particle noise. In the simplest approach the close- 
packed cluster and the discrete-particle noise can be characterized, for each N, 
respectively, by a point-like mass M and mobility Xj and by To, the temperature 
of the circular wall. A stochastic equivalent description is provided in terms 
of the Fokker-Planck equation for P(rn^,t) [V, In the limit F — > 0, the 
(deterministic) steady-state solution of Eq. (|7.17p is Tyji — — ^- the macro- 
particle at rest, located at the center of the box. At F > 0, the steady state 
distribution function is given by the steady state solution of the Fokker-Planck 
equation i2, which is Pfrm.) — (Ttcr^)^^ exp(— r^/ff^), where cr^ = F/k, and 
the normalization constant is computed for ff ^ R. The variance cr^ is the ratio 
of F (a characteristic of the discrete noise) and k (a macroscopic quantity). An 
important insight can be achieved from the N -dependence of cr, obtained by 
molecular dynamics (see Section [7. 2 1) . In analogy with equilibrium systems, one 
might expect the relative magnitude of fiuctuations to decrease with increasing 
N. Surprisingly, this is not what we observed, see Fig. 17.61 One can see that 
ff(N) approaches a plateau, that is fluctuations persist at large N. The small-N 
behavior (the first three data points: N = 150, 315 and 480) in Fig. 17.61 agrees 
with the dependence cr/R = 0(N^^/^), expected for an ideal gas in equilibrium. 
Not surprisingly, for these relatively small N the clustering effect is small. On the 
other hand, at very large N, when the whole system approaches close packing 
(in our event-driven simulations, this corresponds to N = 36275), cr/R must 
go to zero. We could not probe this regime, however, as simulations became 
prohibitively long already at N ~ 15000. An apparently related phenomenon 
has been recently reported for rectangular geometry, in the context of a van der 
Waals-like phase separation in a granular gas driven by a thermal wall |147] . 

^The systematic restoring force which confines the macro-particle has a potential U(rm). 
Not far from the center, U(rm) can be approximated by a harmonic potential. 
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Figure 7.6: Scaled standard deviation of the center of mass position as a function 
of N, obtained in molecular dynamics simulations (dots). The dashed line shows 
is for (j/R = 0.75 N^^/'^. This power law characterizes the fluctuations of an 
ideal (elastic) gas in equilibrium. 



7.5 Conclusions 

We employed a simple model system |155j to investigate the structure and fluc- 
tuations of dense clusters emerging in granular gases. The density profiles, 
obtained in event-driven molecular dynamics simulations, are well described 
by hydrostatic equations which employ the constitutive relations suggested by 
Grossman et al. |156) . This agrees with similar results in rectangular geometry, 
with and without gravity [1301 11561 1161| . We also extended their validity for 
lower restitution coefficient, and for the circular geometry. The agreement is 
particularly good for the large density region (close packing). The azimuthally 
symmetric state is predicted to be stable, fact which is apparent with our sim- 
ulations. The observed Gaussian fluctuations of the center of mass, which are 
too large to be explained by diffusion, suggest an effective Langevin description 
in terms of a macro-particle in a confining potential of hydrodynamic nature, 
driven by discrete particle noise. Surprisingly, the fluctuations persist as the 
number of particles in the system is increased. This unexpected behavior is 
stark contrast with the shown by elastic gases in thermal equilibrium. 

Other important point is the transition from a "smooth state" (with en- 
hanced density in the middle) to a " cluster state" with a pronounced and sharp- 
edged cluster in the center, as a function of particle number density. It has not 
been stated before, that the process of cluster formation requires a certain par- 
ticle number and the transition between a homogeneous (or irregular) state to 
a regularly clustered state (with a persisting cluster in the middle) is rather 
sharp. In our case it is at N = 1700 ± 300. 
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It is worth mention that the global constitutive relations of Grossman et 
al. |156j completely ignore the issue of coexistence of different phases of the 
granulate: the liquid-like phase, the random close-packed phase, etc. However, 
the vibrofluidized steady state, considered in this chapter, has a zero mean flow. 
Therefore the viscosity terms in the hydrodynamic equations vanish. This fact 
is not merely a technical simplification. The shear viscosity of granular flow is 
finite in the liquid-like phase, and infinite in the (multiple) domains of the ran- 
dom close-packed phase. The effective total viscosity of the system is expected 
to diverge when the coarse-grained density slightly exceeds the freezing den- 
sity |145[ 1146] . This invalidates any hydrodynamic description for sufficiently 
dense flows, and requires the introduction of an order parameter and a different 
type of the stress-strain relation into the theory (cf. Ref. |155j ). Luckily, these 
complications do not appear for a zero-mean-flow state. Indeed, the equation 
of state, heat conductivity, and the inelastic heat loss rate do no exhibit any 
singularity around the freezing point, and all the way to the hexagonal close 
packing. Therefore, the hydrostatic description remains reasonably accurate far 
beyond the freezing point. 

Future work should focus on dense clusters (more precisely, on their "core" 
regions) and characterize the positional and orientational ordering at different 
N. This can be achieved by measuring two correlation function: the positional 
correlation function and the orientational correlation function. The motivation 
is the following. There are classical papers of Kosterlitz and Thouless, Young 
and Nelson, and Halperin, from the 70-ies — these yielded the KTHNY theory of 
melting, reviewed by Strandburg in Ref. }119j — , which suggested that melting 
in assemblies of elastic hard disks occurs (as the density decreases) via two 
continuous phase transition: the first one from the perfect hexagonally packed 
phase to the so called /iea;ati(j£| phase, and the second one from the hexatic 
phase to a liquid-like phase. A major question which remains to be answered 
concerning our inelastic system is whether there is any signature of additional 
instabilities, as the density in the central region increases. 



°This phase can exist between solid and liquid phases. The system in the hexatic phase 
has no long-range translational order but has quasi-long-range orientational order. 
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Appendix A 



Equilibrium properties of 
the Lattice Gas 



The driven lattice gas studied in Chapter [3] reduces to the (equihbrium) lattice 
gas [167| for zero field. This is a model for density fluctuations and liquid-gas 
phase transformations. As its nonequilibrium counterpart, this is defined in a 
square lattice with N — L||Lj_ sitefl Thus, each lattice site i can exist in two 
states, occupied by a particle or empty, labelled by an occupation variable rii, = 1 
or rii = 0, respectively. The lattice gas can be mapped to the Ising model with 
ferromagnetic couplings 0. The Ising model is perhaps the simplest system that 
undergoes a nontrivial phase transition in two or more dimensions. The model 
was initially proposed to describe ferromagnetism — the presence of spontaneous 
magnetization in metals such as Fe and Ni below a critical temperature Tc — . 
The relation between the lattice gas and the Ising model is set by simple change 
of variables, namely, Si = (2ni — l)/2, where Si are spin variables capable of 
two orientations, "up" and "down". The spin up Si = 1/2 and the spin down 
Si — —1/2 states correspond to occupied (nt = 1) and unoccupied (nt — 0) 
cells, respectively. Each spin i interacts with one another via an exchange 
interaction Jij > 0. Assuming ferromagnetic coupling and isotropic interactions 
(Jij = J > 0), the Ising Hamiltonian is 



where (ij) denotes neighbor interactions and B is the external magnetic field, 
which plays the role of the chemical potential of the lattice gas. In fact, the 
canonical partition function for the ferromagnetic model is the grandcanonical 
partition function [3] for the lattice gas model. Therefore the lattice model is 

^We consider here only two-dimensional lattices. 

■^Historically, the Ising model 52 was introduced before the lattice model |167| . Likewise, 
the lattice gas as and the binary alloy model 1321 are equivalent. 
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Figure A.l: Phase transition in the Ising model with NN (circles) and NNN 
couplings (squares). Left panel: Spontaneous magnetization M as a function 
of the temperature. The magnetization vanishes abruptly at Tc/Tons = 1 for 
NN couplings, and at Tc/Tons — 2.32 for the NNN case. Otherwise is zero for 
T > Tc in both cases. Right panel: As expected in equilibrium systems [34] the 
order parameter critical exponent do not depend on the dynamics details. The 
solid lines represent slopes of (3ising = 1/8- 

isomorphic with the Ising model. Thermodynamically the descriptions of the 
ising model and the lattice gas are equivalent, but there are simplifying features 
for the Ising magnet. Henceforth we shall refer to them without distinction. 

Lars Onsager [211 [32 gave in 1944 the (exact) partition function for the 
two-dimensional model, in the absence of external magnetic field. The On- 
sager solution showed the existence of a second order phase transition at — 

2/ln ^1 + ^/2j Jkg^ , referred to as the Onsager temperature Tons- The spon- 
taneous magnetization 



is the order parameter for this transition. Above the critical temperature the 
(average) spontaneous magnetization is zero, whereas below criticality M is 
non-zero, thus one finds the system in a ferromagnetic state. Due to the up- 
down symmetry of the Ising system, the spontaneous magnetization is doubly 
degenerate, with an "up" phase and a "down" phase. In the terminology of the 
lattice gas, there should be a dense (liquid) phase coexisting with a dilute (gas) 
phase below the critical temperature. Otherwise, typical configurations are dis- 
ordered. The left panel in Fig. lA.ll shows the second order phase transition for 
nearest-neighbor (NN) couplings and next-nearest-neighbor (NNN) couplings, 
as schematized in Fig. 13.41 in Chapter [31 For NN couplings, the first sum in 
Eq. (lA.ll) runs over the four NN sites, and for the NNN case it runs over the 
eight NNN sites. As is shown in Fig. lA.ll the phase transition occurs exactly 
at Tons for NN couplings. The immediate effect of considering additional neigh- 
bors, as in the NNN case, is the rising of the critical temperature. For NNN 




(A.2) 
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Figure A. 2: Surface plots of the two-point correlation (upper row) and the 
structure factor (lower row). These are for the Ising magnet, or equivalcntly 
the lattice gas, on a Ly x Lj^ = 128 x 128 lattice with NN (left column) and 
NNN interactions (right column) . To avoid complications due to inhomogeneous 
ordered phases, we choose temperatures above criticality. Temperatures are 
T/To„s = 1.09 for the NN and T/To„s = 2.61 for the NNN case. Notice the 
change of color code with each plot. 

couplings, the critical temperature raises to Tc — 2.32 Tons ,94 . In the neigh- 
borhood of the critical temperature the spontaneous magnetization behaves like 
M « (1 — T/Tc)'^^°'°* as T — > Tc — . Here, piling is the order parameter critical 
exponent. The right panel of Fig. lA.ll shows the behavior of M near Tc for NN 
and NNN couplings. In both cases, piling — 1/8. This result is consistent with 
(equilibrium) renormalization group ideas j34j, which set up that the critical 
exponents are universal. It means that there are only a few parameters, such 
as dimensionality and underlying symmetries which are relevant near the criti- 
cal point, therefore critical exponents do not depend on the particular system, 
either magnetic or fluid. 

In zero magnetic field the pair-correlation function is of particular interest 
since it in a sense measures the "degree of order" of the lattice. It is defined in 
general from the partition function as C = (siSj), where (• • • ) (see also Eq. p.6|) 
in Chapter [3]) denotes the proper ensemble average. Two-point correlations 
are shown in Fig. IA.2I Correlations in the equilibrium Ising model are short 
ranged, except at the critical point where the free energy shows a nonanalytic 
point, controlled by a correlation length £,. Specifically, they decay exponen- 
tially (in Fig. |X2] or more clearly in Fig. IA.3I) when r, a typical inter-particle 
separation, becomes large compared to £,. Translated into the momentum space, 
the pair correlation function lead to the structure factor S(Kx, Ky ) (see Eq. (|3.7|) 
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Figure A. 3: Main graph: Parallel C|| and transverse Cj^ projections of the two- 
point correlation function. Symbols o (parallel) and + (transverse) are for NN 
interactions, whereas □ (parallel) and x (transverse) correspond to the NNN 
case. Inset: Parallel Sy and transverse S± projections of the structure factor. 
Here, Tix,ij = 128Kx,ij/27t are integers. Symbols are as in the main graph. 
System size and temperatures are as in Fig. IA.21 



in Chapter [3]), i.e., structure factor is the Fourier transform of the two-point 
correlation function. Structure factors are shown in Figs. IA.2I and IA.3I Expo- 
nential decay in correlations, takes the form of "analyticity at the origin" in the 
momentum space. Notice also that this analyticity in the origin do not depend 
on whether one takes isotropic or anisotropic couplings. As indicated above, we 
have assumed isotropic couplings between spins, and therefore both C and S are 
isotropic (more clearly in Fig. lA.3| ). In Chapter|3]we show that the driven lattice 
gas behaves quite differently, displaying generic power law decay in two-point 
correlations, not only at the critical point, and discontinuity singularity at the 
origin in the structure factor. 



Appendix B 



Interfacial stability of 
Driven Lattice Gases under 
saturating field 



In Chapter [3] we pointed out one of the most striking, even counterintuitive, 
features of the Driven Lattice Gas (DLG) with nearest-neighbor (NN) inter- 
actions: Monte Carlo (MC) simulations with Metropolis rates show that the 
critical temperature Te monotonically increases with the external driving field 
from the Onsager value Tq = Tons = 2.269 Jkg^ to Too — l-4Tons- As we con- 
cluded, this is associated with the fact that a driven particle is geometrically 
restrained in the DLG. In fact, allowing hops and interactions to the next- 
nearest-neighbors (NNN) one observes just the opposite behavior: the critical 
temperature decreases as the field increases. Moreover, there is no phase tran- 
sition for a large enough field, i.e., Te — > as E — > oo. Or in other words, in 
the DLG with NNN interactions a strong field prevents transition to an ordered 
low-temperature phase. This is better illustrated in Fig. l3.2l in Chapter|3l Qual- 
itatively, this difference of behavior in the infinite field limit between the two 
models may be understood in terms of stability of the liquid-gas interface. That 
is, in contrast with the NN (well-ordered) interface cannot be stable 

under a strong field in the NNN case. This suggest a balance between thermal 
effects, which favor the stability (the order) of the interface |101j . and driving 
effects, which would drive particles out of the interface favoring disorder. 

In order to study this balance, the disparate behavior between both dy- 
namics, and characterize the two competing tendencies — i.e, the field vs. tem- 
perature balance — we propose a generalized driven lattice gas model. This is 
an intermediate model between the two cases above mentioned (the NN and 
the NNN case) which captures the behavior in the large field limit observed in 
each system. It is defined on a two-dimensional lattice L|| x Lj^ with periodic 
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boundary conditions in which particles interact via the usual Ising Hamiltonian, 

H = -4^ffj(Jk. (B.l) 

Each site 1 G is either occupied by a particle or empty, which we denote 
by a lattice occupation number Ui taking two values: di = 1 (occupied) or 
(Ji = (empty). The sum in the Hamiltonian runs over all the NNN sites. As 
in the standard DLG, dynamics is induced by the competion between a heat 
bath at temperature T and an external driving field E, which break the detailed 
balance condition. The field is assumed pointing along one of the principal lattice 
directions, say horizontal. Time evolution is then by microscopic dynamics 
according to the Metropolis transition probability per unit time w (this was 
previously defined in Eq. p.2p in Chapter [3]): 

w((j^ ff') =min{l,e(^"+^-^)/'^«T| _ (g_2) 

Here AH is the energy difference of a particle-hole exchange and E • x denotes 
the dot product between the field and the displacement of the exchange, ff = 
{cJi, ; i €E two-dimensional lattice} stands for a configuration. 

Given a particle, its neighboring sites along the principal lattice directions 
has different probability to be reached than the sites along the diagonals (see 
Fig. IB.l| ). The former are always accessible (probability 1), whereas the latter 
have probability < a < 1 , parameter which is set at the beginning of the 
simulation (as well as T and E). The algorithm proceed as follow^: 

1. We randomly select a NNN pair i,j of sites with different occupancies, 
(Ji ^ ffj. 

2. As Fig. IB. II illustrates, if the selected pair is along one of the principal 
lattice directions the exchange attempt is allowed (Go to next item). In 
other case, draw a random number Ti , and if Ti < a then the attempt is 
allowed (Go to next item), otherwise rejected (Go to item 1). 

3. Perform a standard MC trial move controlled by the Metropolis rate 
function w(z) = min{1,e^^}, i.e., controlled by the standard biased rate 
Eq. (jB.2p . To be specific, first compute z; if z < the attempt (particle- 
hole exchange) is finally accepted; if, however, x > 0, we draw a random 
number r2 and perform the exchange only if ti < e^". 

4. Go to item 1. 

Notice that when a ^ only NN particle-hole exchanges are allowed, otherwise 
there is a non-zero probability for diagonal hops. The limit case of a = 1 corre- 
sponds with the NNN dynamics, i.e., with the NDLG model already introduced 

^The algorithm as we use it is more complicated to obtain better (and faster) results. We 
do not present this algorithm due to its complexity. Instead, we show a simpler algorithm 
which would work like this. 
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Figure B.l: Schematic diagram of the accessible sites a particle (at the center, 
marked with a dot) has for the aDLG. Attempts (particle-hole exchanges) along 
the principal lattice directions are always allowed (marked with probability 1). 
Diagonal attempts are allowed with probability a. The limit of a — > 0(1) 
corresponds with NN (NNN) couplings. Once a is set at the beginning, the 
allowed exchanges are controlled by the Metropolis rate function Eq. (jB.2[) . 

in Chapter [3l Therefore, this model (henceforth referred to as the aDLG) can 
capture both the DLG and NDLG phenomenology by varying the parameter oc 
(the effective connectivity). One may also think on a as a temperature which 
controls the diagonal degrees of freedom. For E = this model reduces to the 
equilibrium lattice ga^ (see Appendix IX| . 

We carried out extensive MC simulations on the aDLG by assuming half- 
filled lattices and periodic boundary conditions (toroidal conditions). Hence, 
nontrivial nonequilibrium steady state is set in asymptotically. As in the stan- 
dard DLG, MC simulations reveal that this model undergoes a temperature- 
induced second-order phase transition. This transition is analogous to the stud- 
ied in Chapters |3] and SI However, in this appendix we are rather interested 
in the stability of the interface under saturating fields. Specifically, we study 
the mechanisms which lead to an ordered or disordered low-temperature states 
in the large field limit, i.e., to a stable or unstable interface. To this end, we 
restrict ourselves to the E — > cxd limit and to T — > 0. The latter limit is the 
most favorable case for ordering, which will enable us to identify the point in 
where the field effects outweigh the thermal effects. In addition, this limit is 
convenient for theoretical and computational purposes since T is no longer a 
parameter. 

For a = 0, we observe fully frozen striped configurations, which extends 
along the field direction. Here, a one-strip liquid-like (rich-particle) phase per- 
colates along the field direction (the density of the vapor phase is zero). When 
increasing a, i.e. increasing the weight of the diagonal dynamics, some particles 
are driven out the interface. Nevertheless, if a is not too large, the thermal 
effects still dominate. In such a case, a liquid-like (high-density) phase which is 
striped then coexists with its gas (low density phase). With a further increase 
in a, the thermal effects are overweighted by the field effects and the interface 

^Although, of course, the equilibrium critical temperature depends on a, but not the critical 
properties. The aDLG belongs to the Ising universality class for all a. 
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Figure B.2: The a-dependence of the order parameter m defined in Eq. (jB.3p 
for half-fihed lattices. The main graph shows the results from MC simulations 
for L|| X Lj^ = 256 x 256. The inset shows a zooming of the transition region for 
system sizes L|| x Lj^ = 128 x 128 (dashed line) and L|| x Lj^ = 256 x 256 (sohd 
line). 



becomes unstable. This results in fully disordered configurations (single gas-like 
phase). The transition between ordered and disordered states are rather sharp, 
which suggest a phase transition. This is confirmed. We found a topological 
(a-induced) extraordinary phase transition, which is found to be continuous 
(second order). In order to understand this second order phase transition, we 
introduce an order parameter. One quantity that has been used as a measure 
of order in the two-dimensional DLG and some related systems [IT is 



1 



m : 



2v/p(1 -P) 



(M||)-(Mi) 



(B.3) 



where M|| and Mj^ are squared longitudinal and transverse magnetizations given 

by, 



Ml, 



y 



T 2 



(B.4) 



respectively, where x and y denote the principal lattice directions, and p the 
systems density p = N/(L|[Lj_). The order parameter m is a measure of the 
difference between the densities of the liquid (strip-like) and gas phases: one 
has (M|[) = (M_l) when configurations are fully disordered (as if T = oo), which 
implies m = 0, while (Mj|) — > 1 and (M._l) — > in the limit of T = leading 
to m — > 1 . We computed the order parameter m for different values of a. 
Figure IB. 21 shows how m varies between these two limits as a function of a. As 
observed in Fig. IB. 21 the transition from a stable interface (high ra regime) into 
a unstable one (low m regime) is rather sharp. As a matter of fact, due to the 
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Figure B.3: The main graph shows plots of ra^^'^ for different values of |3: 1/2 
(squares) and 1/6 (triangles); this gives a the straight line only for |3 ~ 0.28 
(circles). The inset is a Log-Log plot assuming (Xc = 0.550 ± 0.004 which gives 
|3 « 0.28. Lines of slope 1/2 (dotted line) and 1/6 (dashed line) are also shown. 



finite size effects some smoothness is introduced in the critical region (see the 
inset in Fig. IB. 21) . It seems natural to try the power law behavior 

m ~ A I ccc — ccl as a ^ ac (from below) . (B.5) 

With this aim, we try to identify a value of (Xc with the familiar Log-Log plots, 
which yields a linear region (see inset in Fig. lB.3|) . The slope near ttc corresponds 
to |3. Alternatively, one may plot m raised to the power of 1/|3 versus a for 
different trials values of (3, looking for also straight lines. This is also shown 
in Fig. IB. 31 The latter procedure has the advantage that no guess for oCc that 
might introduce further errors is involved. Both methods indicate that the data 
are consistent with Eq. (|B.5p . The estimates obtained from this analysis are 
(Xc = 0.550 ± 0.04 and |3 = 0.28 ± 0.02. 

This connectivity-induced (nonequilibrium) phase transition characterizes 
the balance between thermal and field effects, and describes the stability of the 
interface in the aDLG and related models. That is, for low enough connectivities 
the aDLG is well ordered (the thermal effects dominate, and the interface is 
stable; this is indeed the case of the standard DLG), whereas for connectivities 
above the critical point the system exhibits disordered configurations (thermal 
effects are outweighed by field effects, and therefore the interace is unstable, as 
the NDLG). This second order (continuous) phase transition indicates that the 
temperature-field balance is not as trivial as the microscopic dynamic (by the 
Metropolis rate Eq. (jB.2p ) may suggest, i.e., a mere linear superposition, which 
in particular would had led to a linear behavior in Fig. IB. 21 On the contrary, 
this is a clear expression of complexity originated in the cooperative behavior. 
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Appendix C 



Detailed derivation of a 
mesoscopic equation for the 
Driven Lattice Gas with 
next-nearest-neighbor 
Interactions 



In this appendix, we detail the derivation the (mesoscopic) Langevin-type equa- 
tion for the Driven Lattice Gas with next-nearest-neighbor interactions (NDLG) 
introduced in Chapter [3l We empfoy and extend the method introduced in 
Ref. |80j for the Driven Lattice Gas with nearest-neighbor couplings (DLG). 

The original microscopic model consists of a two-dimensional lattice, whose 
sites have occupation variables crt = 1 or 0, for site i. These variables evolve 
following a particle-hole exchange dynamics. Let us define at each point r G 
I?' a density variable, 4)r S which is the averaged value of the occupation 
variables cJi in a region of volume v around r. The system evolves from a given 
configuration y to another y' by choosing at random a particle at point r and 
exchanging it with next nearest neighbor in the \i. direction at time t, namely 



When V is large enough, 4)r is assumed to be a continuous function of r, say 
ct)(r, t). This represent the coarse grained excess particle density field. So that 
we have at time t, 



Cjjr = Cjjr + V ^ (6r - 6r+|x) • 



(C.l) 



y' — {4)(r) + V ^ V^6(r' — r], where cf) e y} . 



(C.2) 



We can also generalized this dynamics to consider exchanges of magnitude y/v 
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with probabiUty ampUtude, the latter being an even function of t) |168| . e.g., 

f(Tl) = (6(i1 + 1) + 6(Tl-1))/2. (C.3) 

Or in other words, Eq. (jC.3P represents the probability distribution function of 
a T], amount of mass attempted to be displaced. The configurations t) evolves 
according to a stochastic hopping dynamics which conserves the number of 
particles, or equivalently the density. That is, configurations have a statistical 
weight P(y, t), which evolves accordingly to the following continuous master 
equation [1] [5]: 



9tP(y,t] =^ 



dTlf(Tl) 



dr [Cl[y Y')P(y, t) - a(y' y)P(y', t)] , (C.4) 



where 0.(y — > y') stands for the transition probability per unit time (transition 
rate) from y to y'. As usual, the transition rates are given in terms of entropic 
and energetic contributions as 

a(y ^ V) = D(AS(ti)) • D(AH(ti) + H,) . (C.5) 

Here D is a function satisfying the detailed balance constraint, i.e., D(— x) — 
e^Dfx), which ensures that in the zero field limit (e = 0) the stationary dis- 
tribution is the equilibrium one, Pst oc e^^. AH and AS are the increment of 
energy H{y') — H(y) and entropy S(y') — S(y), respectively. At a mesoscopic 
level, the equivalent to the microscopic Hamiltonian (Eq. (|3.ip in Chapter [3]) 
is the standard (f)^ (Ginzburg-Landau) Hamiltonian. The structure of the free 
energy consist of two contributions [33] : entropic and energetic which are given 
by, respectively 

S(y)=v|dr (14)^ + ^4)4) 

H(y] =v|dr Q(V4))V^) , (C.6) 

where a and g are entropic coefficients whereas the parameter t comes from 
the energetic functional. Notice that, in Eq. (|C.5|) the increment of energy 
from the drive He enters the dynamics though D(AH(r|) + He), where H^ = 
r| |x • £ (1 — 4)^) + 0[v^^ ). This means that the transition rates depend on the 
energy and entropy difference between configurations plus a term He whose 
dominant part in v^^ is the natural choice to mirror the effects of the drive 
as far as it accounts for the local increment of energy due to the driving field. 
Notice also that the detailed dependence of the coarse-grained field £ on the 
microscopic field introduced in Eq. (|3.2I) remain still as an open issue [11]. We 
also assumed the driving field acting along one of the principal lattice directions, 
say horizontal. 
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The next step is to get a Fokker-Planck equation by expanding the master 
equation Eq. (jC.4p in v^^ up to order (Kramers-Moyal expansion [I]). To 
this end, we expand a functional F around y which read, 



k=l 



•64) 



(C.7) 



where F is a functional representing either S, W, or P. The operator means 
the functional derivative of F(ct)) with respect to cj). To be specific, we also 
choose Hf as 



'|k„k+l.,-k 

H.(y^y)^<-' + I. ^^\V^,^: (^v,^) AU), (C.8) 



with aI/^' = |i,-£ (1 - ct)(r)2). These expansions yield the following Fokker-Planck 
equation [80] , 



9tPt = Y. 



dr ( V,- I X 



(C.9) 



where 



PtZl,Z2j 
q(Zl,Z2) 



Pt = P(y,t), 

dri rif(ri) D(rizi ) D(riZ2), 

dTlTl2f(Tl)D(TlZl)D(riZ2). 



(C.IO) 



After using standard techniques in the theory of stochastic processes [T] we 
derive from the Fokker-Planck equation its stochastically equivalent Langevin 
equation using the Ito prescription, which reads 



at4)(r,t) = [p(A^,A;/ + A^) + q(AS^,A;^ + A^)i/2 f,^(r,t)j . (C.ll) 



The time has been rescaled by v ^ and £,^(r, t] is a delta-correlated Gaussian 
white noise, i.e., (f,^,(r,t))t = and (£,^(r,t) f,^'(r',t')) = 6^^,6(t-t')6(r-r']. 

We focus on the critical region where large fluctuations on all length scales 
dominate. Further simplification in Eq. (jC.lip is possible in this regime by 
dropping the irrelevant terms in the renormalization group sense. Following the 
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standard field theoretic methods let us introduce an external momentum scale 
T and make the following anisotropic scale transformations 



(C.12) 

I'll ^ I'll ) 

Next, we expand the Langevin equation in powers of t around t = 0, kcap- 
ing only the leading terms. For next-nearest-neighbor interactions, the sum Y 
involves a sum for each possible direction, namely Y./ Y.\ which are re- 

spectively the parallel, transverse, and the two diagonals directions. After some 
cumbersome algebra, one can split the terms depending on both diagonal (/ 
and \) into the two parallel and transversal to the driving field £ components. 
The time scale, the transverse noise, and the transverse spatial interaction are 
forced to remain invariant under transformation. With this, one gets z — A and 
5 — (s + d — 3)/2. Different scenarios are now possible depending on the value of 
s . Demanding that the most relevant terms in the parallel and transverse to the 
field scale in the same way, as in the standard analysis of the critical behavior 
of the NDLG (see subsection I3.3.2|) . would lead to s = 2. According to this, 
the Langevin equation for finite (coarse-grained) driving field contains a large 
bunch of £-dependent terms. Similarly to what occurs in Ref. |80| . just taking 
a large enough value of the driving field most of those terms vanish. By doing 
this, one is led to the following Langevin equation 

9t4)(r,t) = -Vl4) + (t+ ^a)Vi4) + + aVfct) + £,(r,t) (C.13) 

The resulting Langevin equation is analogous to the one for nearest-neighbor in- 
teractions (see Eq. (13.101) in Chapter [3]), although Eq. (jC.13l) contains additional 
entropic, not energetic, contributions due to the presence of larger number of 
bonds than in the NN case. 



Appendix D 

Hydrodynamic description 
for a three-dimensional 
granular gas: Clustering 
and symmetry breaking 



A straightforward extension to three dimensions of the circular system studied 
in Chapter [7] is a cylinder in which the curved wall is maintained at a constant 
temperature and the other two are elastic (see Fig. ID.ll) . Particles of diameter d 
move inside the cylinder underging inelastic collisions with a constant restitution 
coefhcient [i. In this appendix we study the clustering and symmetry-breaking 
instabilities employing granular hydrodynamics in such a geometry. To this 
end, we generalize the Grossman et al. |156j constitutive relations to three- 
dimensional situations. 

Constitutive Relations. Before going any further, we derive heuristic consti- 
tutive relations for inelastic gases in three dimensions by employing free volume 
arguments in the vicinity of the close packing, suggesting an interpolation be- 
tween the close-packing limit and the usual dilute-limit relations. Here, we 
follow the strategy taken by Grossman et al. jl56j in two dimensions. 

Let us consider our three-dimensional system, in which there is a temper- 
ature gradient in the radial direction because of the thermalized wall (set at 
temperature Tq). As a result, there will be an energy flux along this direction. 
Consequently — assuming both azimuthal and longitudinal symmetry, i.e, the 
state do not depend neither 9 nor x — , the energy balance equation, which is 
derived from the hydrostatic version of Eqs. (|5.1|) . reduces to 




(D.l) 
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2R 





Figure D.l: A perspective sketch of a cylinder of volume ttHR'^. The curved 
surface is thermalized to temperature To whereas the remaining two surfaces 
are considered elastic. Definition of the cylindrical coordinates, which we will 
employ along this appendix, are also indicated. 



where I is the sink term (mean energy lost per unit volume per unit time) and 
K is the coefficient of thermal diffusivity. Constitutive relations for I and k are 
estimated heuristically following the approach devised in Ref [156] , yielding the 
same expressions as in Eq. ()7.3p . These are 



, and 



I = ^(^-^2^T^T^/^ 

yl 



(D.2) 



where a and y are numerical factors of the order unity. The effect of the system 
dimensionality enters only in a and y, which will differ from the two-dimensional 
case. Concerning the equation of state in the low density limit, as is well-known, 
the ideal gas law holds, P — uT. Contrary, in the high density limit, the mean 
free path I is much less than the particle diameter d, that is, I <C d. Thus, one 
finds 

" • - (D.3) 



Tic (d + l)3 

where ric = Vl/ is the close packing value in three-dimensions |169) . In this 
limit |156) the entropy per particle is S ~ In(l^) + g(T), with g an arbitrary 
function of T. Therefore, employing Eq. (|D.3|) we obtain the pressure in the 
limit n ^ Tie, 

P . ^ (D.4) 
TLc — n 

We therefore propose the following interpolation formula for the pressure 

2n 



P =nT 



Tin -n 



(D.5) 



The mean free path can be also expressed in terms of the density and tempera- 
ture. In the dilute limit one has (in three dimensions) 1 = 1/ {n\/lr\.d?' ) , while 
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from Eq. (ID.3|) one has in the high density hmit 

l«^d. (D.6) 

Again, by using these hmits to interpolate a global expression for the mean free 
path, we find 

; (D.7) 

7tV2nd2 Tic - an 

where a = 1 — ^ . 

The Density Equation. Substituting Eqs. (jD.2p into Eq. (jP.ip . and elimi- 
nating the temperature dependence with Eqs. (|D.5[) and (jD.7l) we arrive to a 
second order differential equation for n: 



r dr 



rJ^ z — 
dr 



= ilS(z). (D. 



This is the density equation. For convenience, the radial coordinate has been 
rescaled by R and a normalized inverse density z(r) — TLc/n is introduced. The 
functions J- and Q read 



Hz) 



[z^ + 4z - 2) [az(z - 1 ) + 27t(z - a)]^ 
(z-a){z-1]V2z3/2(^ + 2]V2 

2^^^^(z-aKz-l)V2 



(z + 2) 



3/27I/2 



and T) = ^^^(1 — |x^) {^)^ ■ The fixed total number of particles N yields a 



normalization condition, which read 

rl 



z-i(r)rdr = f/2, (D.IO) 



where f = ~^2a (l')^ average grain area fraction, and A = H/R the 

system aspect ratio. Equations (jP.Sp and (jP.lOp . together with the boundary 
conditions dz(0)/dr = and z(1] = consi form a complete set. The hydrostatic 
problem is completely determined by three scaled parameters t) , f , and A. Notice 
that the definitions of J^, Q, r\, and f differ from the two-dimensional case (see 
Eqs. (rr6| - (|7lT1) V Figure El shows the numerical solution of Eq. (jDH). This 
is actually a one-dimensional solution, which corresponds with the azimuthally- 
and longitudinally-syimnetTic state. Once the density z(r) = TLc/TL(r) is ob- 
tained, the temperature profile can be determined from the equation of state 
(Eq. (Id3]) ). 



Marginal Stability Analysis. In general, Eq. (|D.8|) provides only one of the 
possible solutions — actually, a one-dimensional solution — . However, as we 
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Figure D.2: Scaled temperature (dashed line) and density (solid line) profiles 
predicted by the hydrostatic theory for r\ — ^0'^ and f = 0.036. 

will demonstrate, this system exhibits symmetry-breaking instabilit}0. That is, 
there are truly two-dimensional solutions in which the longitudinal symmetry 
along X is broken. These solutions can be found by linearizing Eq. (jP.Sp around 
the one-dimensional solution given by Eqs. (jD.8p - (|D.10p — a similar analysis 
was performed in Chapters [5] and [7] — . In general, one can rewrite the energy 
balance equation in terms of the three coordinates r, 0, and x: 

V- (J-(r)Vz) =TiQ(z), (D.ll) 

with T and Q given by Eq. (jD.9p . Substituting z(r, 9,x) ~ z(r) + £cu{r, 9,x), 
assuming uo — S(r)TT(0)O(x), and linearizing Eq. (|D.11|) with respect to the 
small correction e, one obtains: 

(I)(x] = Asin(rax) + B cos(rax) , x e [0, A] 
n(0) = C sin(k0) + D cos(k0] , e [0, 2n) 

r"-hlr'-(^m2 + ^ + ^)r = o, rG[o,i], 

where A,B C, and D are constants, and m and k are, respectively, the longitudi- 
nal and azimuthal wave number determined by the boundary conditions. Here, 
r(r) = E[r)T, Q' denotes the z derivative of Q, and functions and Q are 
evaluated at z = z(r] . For fixed values of rj and f , Eq. (|D.12p and the boundary 
conditions 

O'lO) = O'lA) = (D.13) 

represent a linear eigenvalue problem for k and m. We focus here only in 
the longitudinal wave number m, so that m — j7t/A, where j € {0, N} be- 
cause of the boundary conditions Eq. (jD.13j) . That is, we solve numerically 
the eigenvalue problem for k = 0. If there are nontrivial solutions for the low- 
est mode k = then there will be solutions for higher modes. Moreover, in 



This is in contrast with its two-dimensional counterpart, previously studied in Chapter [T] 
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Figure D.3: Main graph: Marginal stability curve m — m(f) for r\ — 10^. 
The symmetry-breaking instability develops in the parameter region above this 
curve. Inset: Eigenfunction VmSy) corresponding to the eigenvalue ra = 0.585 
forri = 10'^ andf = 2.16-10-^ 

Chapter [7] we proved the stability of the azimuthally-symmetric state — against 
linear perturbations — , hence one expects that the azimuthal wave number k 
does not play a relevant role here. The symmetry breaking instability along 
the longitudinal (not azimuthal) direction is confirmed. Figure ID. 31 shows the 
marginal stability curve m. — Ta(f) found numerically. For a fixed value of r| the 
azimuthally- and longitudinally-symmetric state is unstable (stable) for any m 
below (above) the marginal stability curve. Interestingly, the symmetric state 
remains stable for any m beyond a finite interval of f G (fmin,fmax)- Or in 
other words, for each aspect ratio A > Acrit — Ttj/m., where j = 0, 1 , 2, . . ., the 
symmetric state loses stability with respect the longitudinal mode j |161) . 

Therefore, granular hydrodynamics, employing our heuristic constitutive re- 
lations — an extension of the Grossman et al. approach to three-dimensional 
cases — , may be a good candidate tfor describing the clustering and symmetry- 
breaking instabilities which may occur in inelastic hard-sphere systems. 
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Part IV 

Conclusions 



Conclusions and future 
work 



Unfortunately, despite substantial research, a rigorous statistical-mechan-ical 
description of nonequilibrium instabilities remains elusive. The work presented 
here delves into the nonequilibrium realm, seeking a better understanding of the 
basic features of nonequilibrium phase transitions. Indeed, the study of driven 
diffusive fluids and driven granular gases, besides its technological importance, 
is of general interest because it will contribute to the rationalization of nonequi- 
librium phenomena. The emphasis in our research is on the study — in most 
cases at different levels of description, namely, macroscopic, mesoscopic, and 
microscopic — of the instabilities which occur in those systems. This task was 
approached with both state-of-the-art and novel theoretical as well as compu- 
tational techniques in statistical physics. 

This thesis was structured in four parts. The first Part, comprised of Chap- 
ters m [3l and IH was devoted to the early-time kinetic and steady-state prop- 
erties of driven diffusive fluids; Part (Chapters \5[\6\ and [7]) dealt with both 
the atomistic and macroscopic description of two-dimensional driven inelastic 
gases; Part |llT] provided the appendices (Appendices |X1|B1[C1 and iD)) . Finally, 
Part IIVI summarizes the main results and the original contributions presented 
in this work. We also offer an outlook for principal issues to be investigated in 
future work. 

Chapter [H In the first chapter, we reviewed the general features, phe- 
nomenology, and open issues of the class of systems which we are dealing with 
in the text from a mechanical-statistical perspective. Likewise, we discuss the 
subject of this thesis as well as the main rationale behind it. 

Chapter\^ In this chapter, we presented the model — i.e., the driven lattice 
gas (DLG) — which was the prototype for driven diffusive fluids, and we sum- 
marized some of its already known unusual properties and controversies. The 
relevant background for Chapters |3] and |4] is also provided. This comprises state- 
of-the-art simulation schemes (lattice and off-lattice Metropolis Monte Carlo) 
as well as field-theoretical approaches (statistical field theories). 
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Chapter \3i We have described Monte Carlo simulations and field theoret- 
ical calculations that aim at illustrating how slight modifications of dynamics 
at the microscopic level may influence, even quantitatively, both the resulting 
nonequilibrium steady state and the kinetic properties. With this aim, we took 
as a reference the DLG. We introduced related lattice and off-lattice microscopic 
models in which particles, as in the DLG, interact via a local anisotropic rule. 
The rule induces preferential hopping along one direction, so that a net current 
sets in if allowed by boundary conditions. In particular, we have discussed on 
the similarities and differences between the DLG and its continuous counterpart, 
namely, a Lennard- Jones analogue in which the particles' coordinates vary con- 
tinuously. A comparison between the two models allowed us to discuss some 
exceptional, hardly realistic features of the original discrete system — which has 
been considered a prototype for nonequilibrium anisotropic phase transitions 
(although we claim here is not justified). We outline the main conclusions as 
follows: 

• We found that the continuous, off-lattice model closely resembles the DLG 
in that both depict a particle current, a highly anisotropic liquid-vapor 
interface which extends along the field direction, and the corresponding 
second order phase transition. 

• However, they differ in some essential features. Contrary to the DLG, 
its off-lattice counterpart shows a transition temperature which decreases 
with the field strength. Moreover, unlike for the DLG, there is no phase 
transition for a large enough field. Concerning the early process of kinetic 
ordering, early-time anisotropies in the off-lattice case point along the 
field, contrary to the ones observed in the discrete DLG. 

• These observations suggest that the DLG behavior is highly conditioned by 
the lattice geometry, which acts more efficiently in the DLG as an ordering 
agent than the field itself. Therefore, the striking, counterintuitive feature 
of the DLG that a strong field raises the critical temperature above the 
equilibrium one is a geometrical effect — induced by the lattice topology — 
rather than dynamic. 

• In particular, it ensues that, due to its uniqueness, the DLG does not have 
a simple off-lattice analog. 

• In order to deepen in this situation, we studied the DLG with an infi- 
nite drive extending hops and interaction to next-nearest-neighbors (this 
model was named NDLG). We confirmed that the NDLG behaves closer 
to the off-lattice case: the critical temperature decreases with the field 
and early-time anisotropies point along the field. 

• However, we found that allowing jumping along intermediate directions 
leaves invariant both spatial correlations and criticality. The former was 



117 



monitorized by the two-point correlation function and the equal-time struc- 
ture factor, whereas the latter was determined by standard finite size 
scaling, and corroborated by the anisotropic driven system (ADS) (meso- 
scopic) approach. 

• Therefore, allowing jumping along intermediate directions modifies essen- 
tially the phase diagram — which is not determined by bare symmetry 
arguments — but not features (such as generic power-laws in spatial cor- 
relations) that seem intrinsic of the nonequilibrium nature of the phe- 
nomenon. 

These conclusions are properly complemented by the results in the Chapter |4] on 
criticality in a related off-lattice model. Moreover, we propose that the fact that 
particles are constrained to travel along discrete lattice directions may condition 
not only the properties of the DLG, but also the properties of many other lattice 
models. Further research, under way at present, confirms this on the off-lattice 
versions of the Time Asymmetric Exclusion Process ^W[\57\ (TASEP), and of the 
DLG with nearest neighbor repulsion ^TJ |93] . Further issues include the charac- 
terization of the stability of the interface under strong fields (see Appendix [BJ , 
and the study of the early-time kinetic of phase separation in the NDLG from 
the field-theoretical standpoint, by using the successful ADS approach. 

Chapter Since the DLG was shown in Chapter |3] to be unrealistic in 
some essential sense, we introduced here a novel, realistic off-lattice driven 
fluid, which is a candidate to portray some of the anisotropic behavior in na- 
ture. This is a nonequilibrium Lennard- Jones fluid with an external constant 
driving field, settled on a thermal bath. We studied short-time kinetic and 
steady-state properties of its non-equilibrium phases, namely, solid (mono- and 
polycrystalline) , liquid and gas anisotropic phases. Specifically, we described 
the early-time segregation process as monitored by the excess energy, which 
measures the droplets surface; structural properties of the steady state, namely, 
the radial and azimuthal distribution functions, and the degree of anisotropy; 
transport properties; and also an accurate estimate of the liquid-vapor coexis- 
tence curve and the associated critical indeces. The principal conclusions are: 

• This model seems to contain the necessary essential physics to be useful 
as a prototypical model for anisotropic behavior in nature. 

• This case is more convenient for computational purposes, than others such 
as, for instance, standard molecular-dynamics realizations of driven fluid 
systems. 

• This is the natural extension to nonequilibrium anisotropic cases of the fa- 
miliar Lennard- Jones fluid, which has played an important role in analysing 
equilibrium fluids. Indeed, our model reduces to the (equilibrium) Lennard- 
Jones fluid for zero field. Otherwise, it exhibits a net current, and striped 
structures below a critical point. 
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• Concerning kinetics properties, in spite of the anisotropy of the late-time 
"spinodal decomposition" process, earlier nucleation seems to proceed by 
Smoluchowski coagulation and Ostwald ripening, which are known to ac- 
count for nucleation in equilibrium, isotropic lattice systems and actual 
fluids. 

• Unexpectedly, we have also found that the model critical behavior is con- 
sistent with the Ising, equilibrium one but not with the one for the DLG. 
The main reason for this and other disagreements discussed in this chap- 
ter might be the particle-hole symmetry violation in the driven Lennard- 
Jones fluid. 

• In addition, regarding the modeling of complex systems, one may infer that 
spatial discretization may change significantly not only morphological and 
early-time kinetics properties, but also critical properties. This is in stark 
contrast with the concept of universality in equilibrium systems, where 
critical properties are independent of dynamic details. 

• Finally, we believe this model will motivate new experiments on anisotropic 
spinodal decomposition besides theoretical work. In fact, our observa- 
tions on the early-time separation process and structural properties are 
easily accessible by micro calorimetric and spectroscopy experiments, for 
instance. 

Many other issues remain. Further research, currently in progress, involves de- 
tailed studies of its interfacial properties, the role of the finite size effects, and 
sheared versions of this driven Lennard- Jones fluid. Other questions which we 
would like to raise for future work concern how the particle-hole symmetry and 
the shape of the interaction potential influence the critical properties. A detailed 
analysis of the late-time segregation process, which has already been studied 
both for equilibrium |107l I108j and non-equilibrium cases, including the DLG 
|103|I109) . will be the subject of future investigation. Also useful is the detailed 
study of the structure and morphological features of the (monocrystalline) solid 
and glass-like (polycrystalline solid) phases. The fact that our off-lattice model 
and a related driven-diffusive model [9T| exhibit a critical behavior consistent 
with Ising is rather remarkable, and merits further investigations. 

Chapter O In this chapter, we have provided the framework relevant to 
Chapters [S] and [71 We discussed some unresolved puzzles and recent develop- 
ments in granular fluids, in particular, on the applicability of granular hydrody- 
namics to granular flows. We also presented the basic ingredients of the studied 
models and the computer simulation methods (event-driven molecular dynam- 
ics). 

Chapter We have described hydrodynamic derivations as well as event- 
driven molecular dynamics simulations in a monodisperse granular gas confined 
in an annulus, the inner circle of which represents a "thermal wall" . The quasi- 
elastic limit was considered, and the granular hydrodynamics with Enskog-type 
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transport coefficients [135] was employed. We performed a comprehensive study 
of the clustering, symmetry breaking and phase separation instabilities, which 
enable us to establish a detailed phase diagram. To test and complement our 
theoretical predictions, we performed event-driven molecular dynamics simula- 
tions of this system. This led us to discuss the physics of the instabilities which 
occur in this model. Here we itemize the main conclusions: 

• The zero-flow regime was completely determined by three scaled param- 
eters: the grain area fraction, the inelastic heat loss parameter, and the 
aspect ratio of the annulus. 

• In our event-driven molecular dynamics simulations, azimuthally symmet- 
ric steady states were observed far from the driving wall. In spite of the 
small number of particles employed, these annular states were accurately 
described by the numerical solutions of our granular hydrostatic equations. 

• A marginal stability analysis yielded a region of the (three-dimensional) 
parameter space where the annular state — the basic, azimuthally sym- 
metric steady state of the system — is unstable with respect to small 
perturbations which break the azimuthal symmetry. We computed the 
marginal stability curves and compared them to the borders of the spin- 
odal (negative compressibility) interval of the system. 

• We found that, the physical mechanism of the phase separation insta- 
bility is the negative compressibility of the granular gas in the azimuthal 
direction, caused by the inelastic energy loss. Mathematically, phase sepa- 
ration manifests itself in the existence of additional solutions to the density 
equation in some region of the parameter space. 

• Our event-driven simulations of this system also showed phase separation, 
but it is masked by large spatio-temporal fluctuations. By measuring the 
probability distribution of the amplitude of the fundamental Fourier mode 
of the azimuthal spectrum of the particle density we were able to clearly 
identify the transition to phase separated states in the simulations. 

• We found that the instability region of the parameter space, predicted 
from hydrostatics, is located within the phase separation region observed 
in our simulations. This implies the presence of a binodal (coexistence) 
region, where the annular state is metastable. 

• Clustering, symmetry breaking, and phase separation instabilities should 
be observable in actual experiments. In fact, the model we proposed is 
experimentally accessible, and, therefore, if brought to the attention of 
the appropriate community, might stimulate experimental work on these 
instabilities. 

Of course, it would be interesting to test the predictions of our theory /simulations 
in experiment. By focusing on the annular geometry, we hope to motivate ex- 
perimental studies of these granular instabilities which may be advantageous 
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in this geometry. The annular setting avoids lateral side walls (with an un- 
necessary/unaccounted for energy loss of the particles). A possible experiment 
can employ metallic spheres rolling on a slightly concave smooth surface and 
driven by a rapidly vibrating (slightly eccentric and possibly rough) interior 
circle. While particle rotation and rolling friction may become important, we 
expect that the main predictions of the theory will persist. It would be also in- 
teresting to determine the nature (sub- or supercritical?) of the solutions which 
bifurcate from the azimuthally symmetric state. Further studies may exploit 
the similarities between the phenomenology observed in our model system and 
the one in planetary rings, in which clustering, spontaneous symmetry breaking, 
oscillations, and many other instabilities occur. 

Chapter^ In this chapter we addressed granular hydrodynamics and fluc- 
tuations in a simple two-dimensional granular system under conditions when 
first-principle hydrodynamic descriptions break down because of large density, 
not large inelasticity. We put it into a test — using event-driven molecular 
dynamics simulations — in an extreme case when granular clusters with the 
maximum density close to the hexagonal close packing form (we defined them 
as macro-particles). The model system we considered here, comprises inelasti- 
cally colliding hard disks enclosed by circular boundaries, driven by a thermal 
wall at zero gravity. The main conclusions are summarized as follows: 

• Molecular dynamics simulations showed a sharp-edged close-packed clus- 
ter (macroparticle) of an almost circular shape, weakly fluctuating in space 
and isolated from the thermal wall by a low-density gas. This is due to 
coUisional cooling. The macroparticle formation requires a certain number 
of particles, and, as the particle density is increased, the transition from 
a (quasi-) homogeneous state to a macroparticle is rather sharp. 

• We were able to solve numerically a set of granular hydrostatic equations 
which employ the constitutive relations by Grossman et al. ,156 . We 
found that the density profiles, in a wide range of parameters, agrees al- 
most perfectly — including the close-packed part — with the azimuthally 
symmetric solution of granular hydrostatic equations. We also showed 
that, for the same setting, first-principle Enskog-type relations }135) per- 
form poorly. 

• This agrees with previous results in rectangular geometry, with and with- 
out gravity, and extends the validity of the Grossman et al. relations for 
a slightly lower restitution coefficient, and for the circular geometry. 

• A marginal stability analysis have shown — independently of the consti- 
tutive relations employed — that there are no steady-state solutions with 
broken azimuthal symmetry which would bifurcate from the azimuthally 
symmetric solution. This agrees with the persistence of circular shapes 
for the macroparticle as observed in our simulations. 
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• We addressed, for the first time, fluctuations of the macroparticles, by 
measuring the radial probability distribution function of the center of mass 
of the system. These fluctuations turn out to be Gaussian. 

• The observed Gaussian fluctuations of the center of mass suggest an ef- 
fective Langevin description in terms of a macroparticle in a confining 
potential of hydrodynamic nature, driven by discrete particle noise — i.e., 
the cluster performs a Brownian motion inside the system. The associated 
Langevin equation is provided. 

• In stark contrast with equilibrium systems, in which the relative magni- 
tude of fluctuations decreases with increasing the number of particle^, we 
found that the fluctuations persist as the number of particle in the system 
is increased. 

There are still many open issues. A question that should be addressed is whether 
are there additional instabilities in this model system. For instance, an im- 
portant insight, in analogy with the KTHNY theory |119) of two-dimensional 
melting (for fluids at thermal equilibrium), might be given by the positional and 
orientational ordering of the macroparticles at different number densities. This 
can be achieved by measuring their respective correlation functions. In addition, 
both a general solution of the time-dependent hydrodynamic equations and a 
non-linear stability analysis will contribute to the quest. Other future issue may 
involve the detailed study of instabilities in related three-dimensional systems. 
An outline can be found in Appendix |D] It should be also worthy to study how 
polydisperse granulates behave in this system. 

Despite that the appendices are commonly pushed into the background, we 
consider that our Appendices [B] and [Dl to merit their own conclusions. 

Appendix\B[ In this appendix, we introduced a driven-diffusive lattice model 
(the aDLG) intended to aid in understanding the balance between thermal 
and field effects which occur in these models, such as the DLG and related 
models (specifically the models studied in Chapter [3]). This balance can be 
understood in terms of stability of the interface between the condensed and 
vapor phases (which are elongated along the field direction). Thermal effects 
favor the stability (the order) of the interface, whereas the driving effects would 
drive particles out of the interface, favoring disorder. We studied this dynamic 
balance in the extreme case of the saturating field, zero-temperature limit. The 
main results are: 

• Our model undergoes a second-order phase transition employing as an 
order parameter the measure of the ordering on the lattice as a function 
of the effective connectivity. That is, for low enough connectivities the 
system is well ordered (the thermal effects dominate, and the interface is 

■^For an ideal gas in equilibrium: a ~ ©(N^'^^), where a is the relative magnitude of 
fluctuations and N is the number of particles. 
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stable), whereas for connectivities above the critical point (Xc the system 
exhibits disordered configurations (thermal effects are outweighed by the 
field effects, and therefore the interace is unstable). 

• We computed the critical indexes associated with the transition. 

• This second order (continuous) phase transition indicates that the tempera- 
ture-field balance is not as trivial as the microscopic dynamic may suggest, 
i.e., a mere linear superposition [S3]. On the contrary, this is a clear ex- 
pression of complexity originated in the cooperative behavior. 

Appendix\^ The last appendix addressed clustering and symmetry-breaking 
instabilities in a monodisperse gas of inelastic hard spheres confined in a cylin- 
der, the curved wall of which represents a "thermal wall" . We employed a set 
of granular hydrodynamic equations with heuristic constitutive relations, which 
we derived from free-volume arguments. This model is a simple extension for 
three-dimensional cases of the circular system studied in Chapter [T] 

• Following the Grossman et al. [156\ strategy, we proposed phenomeno- 
logical expressions for the equation of state and the transport coefficients 
for three-dimensional cases. We employed free-volume arguments in the 
vicinity of the close packing, suggesting an interpolation between the 
hexagonal-packing limit and the well-known low-density relations. 

• The stationary solutions of the granular hydrodynamic equations predict 
symmetric clustered states far from the thermal wall. Some of the solu- 
tions bifurcate from the (steady-state) symmetric solution, i.e., hydrody- 
namic also predicts the existence of broken-symmetry states. The pres- 
ence of symmetry-breaking instability in this system is in contrast with 
our findings in its two-dimensional counterpart (see Chapter [7]) in which 
no symmetry-breaking instability was found. 

• Granular hydrodynamics, employing our heuristic constitutive relations, 
may be a good candidate for describing the clustering and symmetry- 
breaking instabilities which may occur in inelastic hard-sphere systems. 
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Resumen en Espanol 



La naturaleza presenta una estructura jerarquica en niveles, con escalas tem- 
porales, espaciales y de energi'a que se extienden desde el mundo microscopico 
al macro scopico. Aunque pueda sorprendernos, a menudo es posible tratar- 
los de forma independiente, e.g., en el caso de gases moleculares. El nivel 
macroscopico, que no es sino el mundo percibido directamente por nuestros 
sentidos, se describe hidrodinamicamente. Esto es, es descrito por funciones 
continuas (o continuas a trozos) de sus coordenadas espaciales r y del tiempo t 
[campos hidrodindmicos) . Por consiguiente, las grandes y exitosas disciplinas de 
la ffsica macroscopica, tales como la mecanica de fluidos, elasticidad, aciistica, 
electromagnetismo, son teorias de campos. Estos campos estan determinados 
por ecuaciones integro-diferenciales que involucran funciones, en principio de- 
sconocidas, de r y t. En cambio, el nivel microscopico describe colecciones de 
un gran niimero de constituyentes, tfpicamente atomos, moleculas o entidades 
mas complejas, los cuales interactiian los unos con los otros de acuerdo con 
ciertas fuerzas electro-mecanicas. Aunque generalmente, la evolucion temporal 
de cada individuo esta dada por las leyes de la mecanica cuantica, en muchas 
situaciones la aproximacion clasica resulta ser excelente. La diferencia entre 
niveles macroscopico y microscopico es esencialmente relativa: el concepto clave 
es el niimero de cuerpos, no su tamancH. 

Si hay una clara separacion entre microscopico y macroscopico, aiin es posi- 
ble introducir una nueva escala intermedia: el nivel de descripcion mesoscopico. 
Esta es una representacion gruesclll que tiene por objeto dar cuenta de la ffsica 
de los modelos microscopicos a escalas temporales y espaciales mas altas, y que, 
matematicamente, esta caracterizada por ecuaciones diferenciales en derivadas 
parciales estocasticas [U [2]. Para ser exactos, la descripcion mesoscopica de 
un sistema dado, cuya dinamica es, en promedio, gobernada por las leyes 
de evolucion temporal macroscopica, es el resultado del recuento y acumu- 
lacion de fluctuaciones microscopicas. Este enfoque nos facilita el estudio de 
las propiedades de variacion lenta de los sistemas de muchos cuerpos y, por 
otra parte, conforma el fundamento teorico de los estudios mas modernos de 
fenomenos cri'ticos. Es comiin encontrar en la literatura trabajos que incluyen 

^Por ejemplo, una galaxia es un objoto (macroscopico) compuesto por un gran niimero de 
estrellas (microscopicas); del mismo modo, una comunidad de seres vivos es considerada como 
un conjunto de individuos; un gas como conjunto de moleculas, etc. 

*0 "coarse-grain" (grano grueso) segiin fuentes bibliograficas en lengua inglesa. 
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en este nivel intermedio descripciones de teoria cinetica de los gases, e.g., la 
ecuacion de Boltzmann. Por ello, este nivel es a veces conocido como cinetic(^ 

Estos tres niveles de descripci6r|f| son, en principio, enfoques equivalentes de 
la misma realidad ffsica. No obstante, las leyes y simetrfas asociadas a cada 
nivel son tan diferentes entre sf que desde muy pronto surgio en la comunidad 
cientifica la necesidad de entender como se relacionan, o, dicho de otro modo, 
como las propiedades hidrodinamicas emergen en terminos de la dinamica mi- 
croscopica de los constituyentes subyacentes. Asi se desarrollo la Mecdnica 
Estadistica, disciplina que proporciona esas relaciones. En efecto, la mecanica 
estadi'stica establece como su propio objetivo el determinar el comportamiento 
macroscopico de la materia originada en el comportamiento colectivo de enti- 
dades (microscopicas) individuales. Algunos de los fenomenos que observamos a 
alto nivel son consecuencia de simples efectos sinergeticos de las acciones de los 
constituyentes, por ejemplo, la presion ejercida por un gas molecular en las pare- 
des del recipiente que lo contiene o miles de luciernagas centelleando al unfsono, 
mientras que otros son ejemplos paradigmaticos de comportamiento colectivo 
emergente, e.g., el cambio de regimen laminar a flujo turbulento en fluidos o el 
millon de atomos que forman el programa de la vida: el acido desoxirribuno- 
cleico o ADN. En estos liltimos casos, el comportamiento de los constituyentes 
llega a ser singular y muy diferente del que tendria en ausencia del resto, dando 
lugar a un comportamiento sin homologo directo en las propiedades o dinamica 
de los constituyentes. 

El mayor logro de la mecanica estadistica es la Teoria de Colectividades [3] . 
Esta establece formalmente el nexo de union entre las propiedades macroscopicas 
de sistemas en equilibrio y las leyes que gobiernan las interacciones a nivel 
microscopico entre partfculas individuales. Se dice que un sistema dado se 
encuentra en equilibrio cuando esta aislado, no presenta histeresis y se haya 
en un estado estacionario en el que todas sus propiedades macroscopicas per- 
manecen fijas [4l. En tal caso, las propiedades macroscopicas son expresadas en 
terminos de variables termodindmicas. En efecto, la termodinamica es, en este 
sentido, una descripcion hidrodinamica que consiste en leyes y relaciones entre 
magnitudes termodinamicas, y cuya base teorica microscopica es la mecanica 
estadistica. Desde el punto de vista matematico, la teoria de colectividades 
se ha axiomatizado por completo. Esto nos proporciona, al menos en princi- 
pio, expresiones analiticas para estas relaciones termodinamicas, permitiendo 
obtener toda la informacion macroscopica relevante del sistema en cuesti6r0. 



^Aunque ecuaciones diferenciales estocasticas y ecuaciones cineticas puedan ser clasificadas 
en un mismo nivel, involucran diferentes escalas (mesoscopicas) tanto espaciales como tempo- 
rales. 

^Nos restringiremos a los niveles aqui descritos, aunque existen escalas hidrodinamicas aiin 
mas altas, e.g., el asi Uamado tren de vortices de Karman. 

^Una vez especificado el Hamiltoniano microscopico del sistema (llamemosle H), la dis- 
tribucion "canonica" de probabilidad de estados estacionaria V esta dada en terminos del 
factor de Boltzmann V = e^"^^^^^ /Z, donde Z es la funcion de particion, kg es la con- 
stante de Boltzmann y T la temperatura. A partir de los adecuados promedios sobre esta, 
los observables estacionarios pueden ser calculados. El resto de dificultades son "meramente" 
tecnicas. 



139 





mmmm 















Figure D.4: La complejidad en la naturaleza se manifiesta a todas las escalas. 
Algunos ejemplos Uamativos son los siguientes. Arriba izquierda: geometrfas 
fractales en cuencas fluviales en El Yemen (fuente NASA). Arriba derecha: mod- 
elado eolico franjiforme en las Great Sand Dunes National Monument, EE.UU. 
(fotografiado por Bob Bauer). Abajo izquierda: tornado como un flujo geofisico 
en rotacion. Abajo derecha: oscilones en suspensiones coloidales agitadas ver- 
ticalmente [S]. 

por ejemplo, la energia libre, la ecuacion de estado, la radiancia espectral, etc. 

Por contra (jy afortunadaniente!), en la naturaleza los sistemas en equilibrio 
son la excepcion (incluso una idealizacion) mas que la regla: fenomenos fuera 
del equilibrio son, con creces, mucho mas abundantes. Galaxias, seres humanos, 
reacciones quimicas, flujos geofisicos, portadores de carga en dispositivos semi- 
conductores, mercados financieros, trafico en autopistas, por citar unos pocos, 
son sistemas de muchos cuerpos en condiciones de no-equilibric@. Los sistemas 
fuera del equilibrio se caracterizan por no estar cerrados, es decir, por intercam- 
biar energia, particulas y/o informacion con su entorno. En general, el estado de 
un sistema fuera del equilibrio no es determinado linicamente por ligaduras ex- 
ternas, sino que tambien depende de su historia (presentan histeresis). Esto da 
lugar a la enorme complejidad del mundo real, la cual surge a todos los niveles 
de descripcion [SJ [71 IHl [TU] y se manifiesta en formacion de patrones, fractales, 
sincronizacion, caos, criticalidad auto-organizada, etc. Toda esta sorprendente 
y complicada fenomenologia es comiinmente asociada a inestabilidade^, las 

* Fuera del equilibrio se encuentra tambien cualquier sistema que se encuentre relajando 
hacia un estado de equilibrio. No obstante, en la mayon'a de estas situaciones su dinamica 
puede comprenderse adecuadamente a partir de conceptos de equilibrio. 

^Las inestabilidades presentes en sistemas lejos del equilibrio presentan similitudes a nivel 
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cuales se describen como cambios de fase [3 [TTJ [121 ES] , bifurcaciones, sinergia, 
etc. con el proposito de conectar lo microscopico con las estructuras coherentes 
observadas a alto nivel. 

Sin embargo, la mayoria de los estudios de sistemas fuera del equilibrio real- 
izados hasta la fecha adoptan un punto de vista eminentemente fenomenologico 
(macroscopico). En general, poco se sabe del porque de tales estructuras o 
de como esta complejidad emerge desde interacciones colectivas a nivel mi- 
croscopico. Consecuentemente, no existe una teoria apropiada y el desarrollo 
de solidos fundamentos matematicos esta solo en sus primeros pasos comparado 
con el del equilibrio, donde la teorfa de colectividades funciona con exito. 

Una diferencia fundamental entre mecanica estadfstica del equilibrio y del 
no-equilibrio es que, mientras en la primera la distribucion de probabilidad de 
estados, de la cual se dcriva toda la informacion macroscopica litil, es conocida, 
fuera del equilibrio es, a priori, desconocida. Uno debe encontrar una dis- 
tribucion de estados, en este caso dependiente del tiempo, la cual vcrifica una 
ecuacion general de evolucion, a saber, la ecuacion maestra [T]. Los casos resol- 
ubles se limitan a unos pocos, y lo son solo de forma aproximada. TJnicamente 
es posible ofrecer un enfoque unificado para sistemas que se encuentran no 
demasiado lejo^^ del equilibrio |14| , donde el sistema puede ser tratado pertur- 
bativamente alrededor del estado del equilibrio empleando las tecnicas clasicas 
de las teorias de respuesta lineal. En cualquier caso, nuestra atencion se cen- 
tra precisamente en sistemas lejos del equilibrio, donde tales enfoques dejan de 
tener utilidad. 



En este contexto, el objeto principal de esta tesis reside en el 
estudio de inestabilidades — a diferentes niveles de descripcion — 
en sistemas forzados, i.e., mantenidos lejos del equilibrio par un 
agente externo. Nos centraremos en das importantes clases de sis- 
temas, a saber, fluidos difusivos con arrastre y gases granulares vibro- 
fluidificados. 

La primera familia (fluidos difusivos con arrastre) incluvj"! sistemas que se 
encuentran acoplados a dos fuentes de energia constante, de tal modo que se 
establece un flujo estacionario de energia a traves de aquellos. Esta definicion 
es ciertamente bastante amplia e incluye sistemas de muy diversa indole. Aqui 
nos restringiremos a sistemas que presentan una corriente de particulas (cuyo 
niimero es una magnitud conservada) a su traves, anisotropias espaciales asoci- 
adas a la accion de un campo externo y en los cuales finalmente se alcanza un 
estado estacionario de no-equilibrio. Un ejemplo sencillo y prototfpico podria 



morfologico con las observadas en el equilibrio, descritas estas de forma mecano-estadistica. 
Por tanto, suelen referirse como cambios de fase del no-equilibrio, aunque los cambios de fase 
fuera del equilibrio, sin las restricciones de este, son mucho mas ricos desde un punto de vista 
fenomenologico. 

^"Como es de suponer, una definicion precisa de "no demasiado lejos" es muy difi'cil de dar 
y depende en gran medida del sistema estudiado y debe darse ad hoc. 
^^Aqui seguimos la definicion dada por Schmittmann & Zia en Ref. |15j . 
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ser el de una impedancia en un circuito electrico disipando a la atmosfera la 
energia suministrada por una bateria. Pero incluso en esta reducida clase de 
sistemas, la distribucion estacionaria de probabilidad es desconocidcF^. 

La motivacion fundamental subyacente tras estudio de fiuidos difusivos con 
arrastre es la necesidad de desentranar los ingredientes basicos del compor- 
tamiento observado en sistemas lejos del equilibrio. Por ello, es justificado, 
pertinente y necesario el estudio de modelos simplificados que capturen lo esen- 
cial del comportamiento microscopico que da lugar a las complicadas estructuras 
macroscopicas. Estos modelos son a menudo caricaturas de los sistemas reales 
en los cuales la dinamica de los constituyentes esta dada por simples reglas es- 
tocasticas. Entre ellos, los modelos reticulares [TTl [TSl HH HZl [THl [TH [5D] han 
jugado un papel preponderante debido primordialmente al hecho de que, a ve- 
ces, aportan resultados anali'ticos exactos y permiten aislar las caracterfsticas 
esenciales de un cierto sistema. Ademas, los modelos reticulares nos ayudan a 
ganar la intuicion necesaria para desarroUar teorias y, desde un punto de vista 
experimental, su implementacion en computadores resulta ser muy eficiente. 
Todo esto ha contribuido al desarrollo y aparicion de numerosas y novedosas 
tecnicas, incluyendo la teoria estadistica de campos lejos del equilibrio. En con- 
creto, modelos reticulares han modelado con exito algunos aspectos esenciales 
de organismos con estructura social [211 ESI El], formacion de redes fluviales 
[21], propagacion de epidemias [53], vidrios, circuitos electricos, enzimas [5S], 
trafico [26l|27l[29], coloides y espumas, mercados financieros [30], etc. 

Los cambios de fase (de equilibrio) y los fenomenos crfticos son actualmente 
bien entendidos, en muchos aspectos, gracias a la aplicacion de los metodos del 
grupo de renormalizacion [3311311 a modelos reticulares. El resultado primordial 
derivado del estudio de fenomenos crfticos (en equilibrio, insistimos) es el con- 
cepto de univers alidad: en las cercanfas de los puntos cn'ticos el comportamiento 
de sistemas dispares esta determinado linicamente por ciertas caracterfsticas co- 
munes basicas — dimensionalidad, el rango de las inter acciones, simetrfas, etc. — 
y no parece depender en gran medida del sistema concreto [33] [351 |3S] . 

La pregunta de si este concepto es tambien aplicable a cambios de fase lejos 
del equilibrio esta todavfa lejos de ser respondida. Por ello, se espera que, al 
igual que para sistemas en equilibrio, los modelos reticulares jueguen tambien un 
importante papel aquf. Trataremos con detalle estos aspectos en el CapftulojSl 
en el que estudiaremos un modelo reticular difusivo con arrastre, paradigma de 
cambios de fase anisotropicos lejos del equilibrio. Sin embargo, es facil imaginar 
que los modelos reticulares, comparados directamente con sistemas reales, pre- 
sentan la desventaja de que son demasiado "crudos" o simplificados. Esto debe 
entenderse en el sentido de que los modelos reticulares no capturan todos los 
detalles, especialmente los morfologicos, del diagrama de fases. Discutiremos 
este y otros problemas en el CapftulolH donde presentaremos un nuevo modelo 
mas realista para simulacion de fiuidos anisotropicos. 

^•^Es necesario hacer notar que en la literatura estos sistemas son referidos como sistemas 
difusivos con arrastre en vez de fiuidos difusivos con arrastre, denominacion que adoptamos 
en esta tesis. Como explicaremos con cierto detalle en el Capitulo [S] reservamos aquella 
denominacion para una descripcion mesoscopica de campo medio. 
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La segunda clase de sistemas considerados en esta tesis versa acerca de gases 
granulares agitados mecdnicamente. Un medio granular [37] es un conglomerado 
de muchas particulas de tamano macroscopico, y por lo tanto, susceptible de ser 
modelados clasicamente. Hay ejemplos por doquier y se presentan a muy difer- 
entes escalas: desde materiales polvoreados hasta nubes de polvo intergalactico. 
Dentro de ese intervale podemos encontrar arena, hormigon, cereales, flujos 
volcanicos, los anillos de Saturno y otros muchos ejemplos. Su importancia 
cientifica y tecnologica es indudable, con numerosas aplicaciones farmaccuticas, 
en la construccion e ingenieria civil, quimicas, agrarias y alimentarias, y, a mayor 
escala, en procesos geologicos y astroffsicos. Ademas, su estudio esta correla- 
cionado con otras clases de sistemas como coloides, espumas y vidrios. Asf pues 
la fenomenologi'a mostrada por estos sistemas es amplisima, yendo desde unas 
pocas micras a miles de kilometres: formacion de patrones en granulados fluid- 
ificados [351 IMl HI] , agregacion [HI , flujos y bloqueos en silos y tolvas [JS] , 
avalanchas y deslizamientos [35] , mezcla y segregacion [37] , conveccion [351 US] ) 
erupciones [SOIET], por citar unos pocos. Por todo ello, su importancia esta 
lejos de toda duda y la adecuada comprension de sus propiedades no es solo una 
necesidad industrial urgente, sino que tambien supone un reto importante para 
la fisica. 

Puesto que los granos individuales tienen dimensiones macroscopicas, la 
friccion, deformacion, ruptura y las perdidas energeticas por colision producen 
disipacion de energfa, esto es, la energi'a cinetica es transferida continuamente al 
medio en forma de calor. Asimismo, la escala de energia tipica para, por ejemplo, 
un grano de arroz de, digamos, longitud I « 5mm y masa m es su energi'a po- 
tencial mgl, donde g es la aceleracion gravitatoria. Es para tamafios de grano 
mucho mas pequefios cuando otros tipos de interaccion ganan en importan- 
cigj^. De ahi que en los medios granulares, y a diferencia de gases moleculares, 
la temperatura T no juega ningiin papel, es decir, keT <C mgl. Los efectos 
dinamicos pesan mas que cualquier tipo de consideraciones entropicas [37j y, 
consecuentemente, argumentos mecano-estadfsticos o termodinamicos resultan 
de poca utilidad. Ademas, dado que cualquier granulado esta formado por un 
mimero finite de constituventej^. frecuentemente el efecte de las fluctuacienes 
en observables macroscopicos es enorme y, por tanto, su estadfstica es dominada 
por los, asf llamados, eventos raros. Otros efectos tambien importantes estan 
relacionados con el fluido intersticial, como aire o agua, aunque en muchas situa- 
ciones y, en buena aproximacion, puede ignorarse considerando linicamente las 
interacciones de grano con grano. En esta tesis (concretamente en la Parte lll| 
nos restringiremos al estudio de sistemas granulares secos, donde la interaccion 
entre particulas predomina sobre la interaccion con el posible fluido intersticial. 

Sin embargo, y a pesar de la dinamica clasica que los caracteriza, los medios 

^''En granos con d ~ 0.1 aparecen cargas electrostaticas superficiales, o para d ~ 0.01 los 
efectos magneticos y de adhesion superficial comienzan a ser relevantes. 

^■^A diferencia de gases moleculares, N <^ donde N es el numero de particulas y N;\ 
el numero de Avogadro. No obstante, si N fuera lo suficientemente grande (a determinar ad 
hoc) aiin seria posible emplear un tratamiento estadistico. 
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granulares se comportan de forma poco convencional: sus fases — solida, Hquida 
y gaseosa — ponen de manifiesto una compleja fenomenologia de caracter colec- 
tivo que los distinguen fuertemente de sus honiologos moleculares. El resultado 
as una pauperrima comprension de sus propiedades dinamicas y estadi'sticas. De 
hecho, desde un punto de vista teorico, solo existen predicciones de caracter muy 
general y muchas cuestiones permanecen sin responder. Los solidos granulares 
son altamente histereticos, muestran indeterminacion de tensiones estaticas, y 
el entramado de granos da lugar a cadenas de fuerzas, tensiones y esfuerzos in- 
ternos, bloqueos, estres conducido, dilatacion o compresion por cizalladura. Los 
liquidos son viscosos, presentan avalanchas, segregacion por tamafios, bandas 
de cizalladura, formacion de patrones. Los gases, los cuales solo perduran con 
un continuo aporte de energi'a (por ejemplo, agitacion mecanica^^l). son alta- 
mente compresibles y muestran agregados inhomogeneos, colapso inelastico en 
modelos computacionales, distribuciones de velocidad no maxwellianas, falta de 
separacion entre escalas macroscopica y microscopica. 

Toda esta fenomenologia hace que resulte muy diffcil enmarcarlos dentro de 
la clasificacion traditional de estados de agregacion de la materia, dificultando 
considerablemente su estudio. En Capitulos[n]y[7]abordamos algunos problemas 
no resueltos y controversias, a saber: cuales son las propiedades estadisticas de 
los medios granulares, como y bajo que condiciones cambian de fase y cuales 
son los modelos continuos optimos y cuando aplicarlos. Especialmente nos cen- 
tramos en la descripcion hidrodinamica de las inestabilidades de agregacion 
heterogenea y ruptura espontanea de simetrfa en gases granulares fluidificados. 

Perspectiva 

A lo largo de esta tesis abordamos el estudio de estas dos familias de sis- 
temas. Nuestro trabajo se divide en cuatro partes. En el Capitulo[l]se presenta 
una introduccion general que describe el objeto principal de esta tesis. Los 
Capitulos m E] y 0] conforman la Parte U que esta dedicada a fluidos difusivos 
con arrastre. La Parte trata de gases granulares y reiine los Capitulos [5l [6] 
y [T] La Parte IIIII esta formada por los apendices El [B] [C] y El Finalmente la 
Parte HVl cierra con un resumen y las principales conclusiones. 

La primera parte comienza en el Capitulo |2l presentando el Gas Reticular 
con Arrastre, conocido como Driven Lattice Gas [111 1151 153] — modelo pro- 
totipico para el estudio de cambios de fase anisotropicos lejos del equilibrio. Las 
partfculas que lo forman se mueven en un reticulo e interaccionan de acuerdo 
con ciertas reglas locales de caracter anisotropico. Estas inducen un flujo de 
particulas en una direccion de tal modo que, si las condiciones de contorno asf 
lo permiten, se establece una corriente. Este modelo es considerado una simpli- 
ficacion de ciertos problemas de trafico y fluidos. Presentamos algunas de sus 
ya conocidas propiedades y tambien controversias, asf como los fundamentos de 
los metodos teoricos y computacionales empleados en los Capitulos El y IH 



Se habla cntonccs de gases fluidificados o vibrofluidificados 
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Una de las preguntas que cabe plantearse sobre cualquier modelo fisico es si 
el comportamiento que manifiesta es universal. En el Capitulo [3] abordamos 
esta importante cuestion estudiando la segregacion anisotropica de fases en el 
modelo prototi'pico introducido en el Capi'tulo [H El enfasis de nuestro estudio 
esta puesto en la influencia de los detalles de la dinamica microscopica en el 
estado estacionario del no-equilibrio resultante, aunque tambien se presta cierta 
atencion a aspectos cineticos. En especial trataremos, tanto desde un punto de 
vista computacional como teorico, las similitudes y diferencias entre el driven 
lattice gas y otros modelos reticulares y no reticulares relacionados con aquel. 
La comparativa nos permite discutir ciertas propiedades excepcionales y poco 
realistas del modelo reticular original. Ademas, ponemos a prueba la validez de 
dos teorias mesoscopicas actuales que pretenden dar cuenta del comportamiento 
cn'tico de esta clase de sistemas. 

En cl Capitulo [4] introducimos un novedoso modelo computacional de 
fluido no reticular y de no-equilibrio disefiado para el estudio de fenomenos 
anisotropicos. Su dinamica e interacciones son mas realistas que en, por ejem- 
plo, el driven lattice gas. Se trata de un modelo muy accesible desde el punto 
de vista computacional, que exhibe una corriente neta de particulas y estruc- 
turas franjiformes a bajas temperaturas orientadas segiin la direccion de un 
campo externo, asemejandose a muchas situaciones reales. Describiremos tanto 
las propiedades cineticas a tiempos cortos (nucleacion) como las estacionarias 
del modelo; en especial nos centraremos en su diagrama de fases: identificamos 
y caracterizanios sus fases solida (mono- y policristalina) , h'quida y gaseosa, to- 
das ellas altamente anisotropicas. El calculo preciso de sus parametros crfticos 
nos lleva a discutir el papel de las simetrias en el modelado de fluidos fuera del 
equilibrio con teon'as de campos y simulaciones. Ademas se esbozan algunas 
comparaciones con experimentos reales. 

El Capitulo [5] sirve de introduccion a la PartelUde esta memoria; parte que 
es dedicada a los gases granulares agitados mecanicamente o vibrofluidificados. 
Perfilamos el estado del arte con respecto a la investigacion en gases granulares. 
Concretaniente discutiremos la aplicabilidad de las teon'as hidrodinamicas al 
estudio de flujos granulares. Tambien proporcionamos brevemente los ingredi- 
entes basicos de los modelos estudiados en los capi'tulos siguientes y los metodos 
de simulacion. 

En el Capitulo [6] cmpleamos una descripcion liidrodinamica granular para 
estudiar la fenomenologfa observada en un gas inelastico confinado en un recipi- 
ente anular sin gravedad. Nuestro principal proposito es caracterizar apropiada- 
mente sus estados estacionarios (de no-equilibrio) y establecer de forma precisa 
un diagrama de fases mediante la descripcion hidrodinamica y la computacional. 
Empleamos relaciones constitutivas de tipo Chapman-Enskog, derivadas estas 
de primeros principios. Prestamos especial atencion a las inestabilidades ob- 
servadas, a saber, separacion de fases, agregacion heterogenea de particulas y 



145 



ruptura de simetria. Esta ultima constituye una prueba precisa para modelos de 
flujos granulares y fomenta el estudio de formacion de patrones lejos del equilib- 
rio. Con este estudio pretendemos motivar experimentos en inestabilidades de 
separacion de fases en gases granulares, que podn'an ser faciles de implementar 
en esta geometn'a. Hasta ahora este tipo de estudios ha empleado geometn'as 
planas en vez de radiales. 

El objeto del Capitulo [7] es el estudio de hidrodinamica granular y fluc- 
tuaciones en un gas granular bidimensional a densidades elevadas, condicion 
esta que impide la aplicabilidad de descripciones hidrodinamicas derivadas de 
primeros principios. Estudiamos un modelo de discos rfgidos cuasielasticos in- 
teractuando en un recipiente circular y ganando energia del contorno. Me- 
diante enfriamiento por colision entre partfculas se forma en la region cen- 
tral un granulado de alta densidad, el cual se comporta como una macro- 
particula. Algunas caracteristicas de esta macropartfcula son bien descritas 
por la solucion numerica de un conjunto de ecuaciones hidrostaticas granu- 
lares, las cuales hacen uso de relaciones constitutivas mas fenomenologfcas. Las 
predicciones hidrostaticas se comparan de forma excelente con simulaciones 
de dindmica molecular por eventos. Ademas desarroUamos una descripcion 
mesoscopica (tipo Langevin) para la macropartfcula, confinada por un potencial 
armonico de naturaleza hidrodinamica y forzada por un ruido bianco de niimero 
finito. 

La Parte inil esta compuesta de cuatro apendices: en el Apendice[A]sc com- 
pendian las propiedades mas relevantes del Gas Reticular (o Lattice C?aj^. el 
cual es el homologo de equilibrio del driven lattice gas tratado en el Capitulo [Sj 
en el Apendice [B] introducimos un modelo cuyo proposito es caracterizar la 
estabilidad de la interfase en modelos reticulares con arrastre; Apendice [C] 
detalla las derivaciones de teorfa estadfstica de campos desarrolladas en el 
Capitulo |3] concernientes a ecuaciones de tipo Langevin; y finalmente en el 
Apendice [D] proponemos un conjunto de relaciones constitutivas para la de- 
scripcion hidrodinamica de gases granulares tridimensionales. 

Finalmente, tras los apendices, se incluye un amplio resumen (que adjun- 
tamos a continuacion) de las principales conclusiones, resaltando las contribu- 
ciones originales de esta tesis y seiialando las Ifneas maestras a seguir en trabajos 
futures. 



Conclusiones 

Desafortunadamente, y a pesar de los muchos esfuerzos cientfficos realizados, 
la descripcion mecano-estadistica rigurosa de inestabilidades lejos del equilibrio 
sigue siendo esquiva. Los trabajos presentados en esta tesis suponen un nuevo 



Modelo de Ising conservado. 
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paso adelante para entender los cambios de fase lejos del equilibrio. En efecto, 
el estudio de fluidos difusivos con arrastre y gases granulares vibrofluidizados, 
aparte de su importancia tecnologica, ahonda en este aspecto y contribuye a la 
racionalizacion de los fenomenos del no-equilibrio. Se ha hecho especial enfasis 
en el estudio y caracterizacion de las inestabilidades — en la mayoria de los 
casos a distintos niveles de descripcion, a saber, macroscopico, mesoscopico y 
microscopico — que tienen lugar en aquellos sistemas. Esta labor ha sido Ue- 
vada a cabo empleando modernas y novedosas tecnicas en el ambito de la ffsica 
estadfstica. 

Capitulol^ Se presenta una introduccion que incluye una descripcion general 
de los sistemas estudiados, su lugar en la ffsica y en la tecnica, fenomenologia 
y problematica. Asimismo se resaltan las motivaciones y el objeto principal de 
esta tesis. 

Capitulo [3 En este capitulo presentamos el modelo considerado hasta la 
fecha prototipo de fluido difusivo con arrastre: el Gas Reticular con Arrastre 
o Driven Lattice Gas (DLG). Repasamos tanto sus propiedades mas relevantes 
como las mas controvertidas. Asimismo, se proporcionaron las bases teoricas y 
computacionales necesarias para los subsiguientes estudios de los Capi'tulos [3] y 
m estas son modernas tecnicas computacionales y teorfas estadi'sticas de campos. 

Capitulo\^ Describimos simulaciones tipo Monte Carlo y desarrollos analiticos 
de teoria de campos que pretenden ilustrar como, a diferencia de lo que ocurre 
en equilibrio, ligeras modificaciones en la dinamica microscopica pueden influen- 
ciar, incluso cuantitativamente, tanto los estados estacionarios como la cinetica 
de relajacion. Con este proposito, se tomo como referenda el DLG introducido 
en el capitulo anterior. Presentamos modelos reticulares y no reticulares rela- 
tives al DLG en los cuales las particulas, al igual que en el DLG, interaccionan 
de acuerdo con ciertas reglas locales anisotropicas. Estas reglas facilitan el 
movimiento preferente en una direccion espacial, de tal modo que se establece 
una corriente a traves del sistema si asi lo permiten las condiciones de con- 
torno. En concrete, se discutio las similitudes y diferencias entre el DLG y su 
"analogo" no-reticular, a saber, un modelo en el cual las coordenadas de las 
particulas varian de forma continua y en el que estas interaccionan entre sf con 
un potencial atractivo de tipo Lennard- Jones 12-6 (tambien conocido como po- 
tencial Weeks-Chandler- Andersen). La detallada comparacion de la morfologia 
observada nos Uevo a considerar ciertas propiedades del DLG como poco realis- 
tas y excepcionales, a pesar de haber sido considerado prototfpico. Concluimos: 

• El modelo continue no-reticular presenta, al igual que su analogo el DLG, 
una corriente neta no nula de particulas, una interfase fluido-vapor alta- 
mente anisotropica (franjiforme) y un cambio de fase de segundo orden 
(continue). 

• Sin embargo, estos difieren en algunas propiedades de caracter esencial. 
Contrariamente a lo que ocurre en el DLG, su analogo no-reticular exhibe 
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una temperatura critica que decrece con la intensidad del campo externo. 
Es mas, a diferencia del DLG, no existe ningiin cambio de fase en el limite 
de campo externo saturante. En lo concerniente a los procesos cineticos 
de formacion de estructuras macroscopicas u ordenacion (nucleacion) , se 
observaron agregaciones de parti'culas con forma triangular apuntando en 
el sentido del campo externo, justamente al contrario que en el modelo 
discreto. 

• Estas observaciones sugieren una falta de universalidad en el modelo retic- 
ular originario. Demostramos asf que el comportamiento del DLG esta 
condicionado en gran medida por la geometrfa del reticulo. Este es el 
agente que favorece la ordenacion, y no el campo en sf mismo. Por tanto, la 
caracten'stica sorprendente y contraintuitiva de que en el DLG "el campo 
externo" eleva la temperatura critica respecto a la del equilibrio no es sino 
un efecto geometrico o topologico y no dinamico. 

• De aquf se deriva, en concreto, que el DLG no tiene un homologo no- 
reticular (mas "realista") inmcdiato. 

• Para profundizar en esta situacion, se estudio el DLG bajo la accion de 
un campo externo saturante a la vez que se aumento el rango de las inter- 
acciones a nuevos vecinos. El modelo resultante, definido como ndlgEzI, 
se asemeja mas al del caso no-reticular: la temperatura critica decrece 
con el campo y las anisotropias triangulares observadas en el proceso de 
nucleacion a tiempos cortos se orientan en el sentido del campo. 

• No obstante, mostramos que, aiin ampliando el radio de interaccion entre 
partfculas, tanto el decaimiento potential de las correlaciones espaciales 
como la criticalidad permanecen invariantes. Las primeras fueron moni- 
torizadas por la funcion de correlacion a dos cuerpos y por el factor de 
estructura, mientras que las propiedades criticas se determinaron emple- 
ando tecnicas estandar de analisis de tamano finito y corroboradas por la 
teoria mesoscopica de campos conocida como anisotropic driven system 
(ADS). 

• Por consiguiente, al extender la dinamica en el DLG se modifica de forma 
esencial el diagrama de fases — el cual no puede ser determinado por crudos 
argumentos de simetrias — pero no las propiedades intrfnsecas relacionadas 
con la naturaleza de no-equilibrio del fenomeno — como el decaimiento de 
las correlaciones. 

La conclusiones de este capitulo estan complementadas apropiadamente por los 
resultados de criticalidad del modelo affn expuesto en el CapitulolH Una posible 
y recomendable continuacion del trabajo desarrollado en este capitulo es el estu- 
dio de la segregacion de fases en el modelo NDLG desde el punto de vista de la 
teorfa de campos empleando el enfoque de la teoria ADS. Otros problcmas son 



DLG con interaccion a vecinos siguientes o next-nearest-neighbors DLG. 
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el estudio detallado de la estabilidad de la interfase bajo campos intensos (vease 
Apendice [B|) y la extension de nuestros metodos a otros modelos reticulares 
como, por ejemplo al Time Asymmetric Exclusion Process (TASEP) [16 ( 157 ) . 

Capitulo Puesto que, como fue demostrado en el capitulo anterior, el 
DLG es poco realista en ciertos aspectos esenciales, presentamos aquf un modelo 
novedoso y mas realista de fluido con arrastre, ideal para el estudio de fenomenos 
anisotropicos. Se trata de un modelo no-reticular lejos del equilibrio en el que 
las partfculas, que interactiian entre si via un potencial tipo Lennard-Jones 12-6, 
son conducidas por la accion de un campo externo, cuya energia se disipa en un 
bafio termico. Estudiamos la cinetica de separacion de fases a tiempos cortos 
y sus diferentes fases estacionarias (altamente anisotropicas), a saber, solido 
(mono- y policristalino), liquido y gas. En concreto, describimos los procesos 
de nucleacion, monitorizados por el exceso de energia, el cual da cuenta de la 
energia interfacial; las propiedades estructurales de las diferentes fases en el 
estado estacionario, caracterizadas por las funciones de distribucion radial y 
azimutal y el grado de anisotropia; y propiedades de transporte. Elaboramos 
tambicn una precisa estimacion de la curva de coexistencia liquido-vapor y de 
los parametros criticos asociados. Las principales conclusiones son: 

• El modelo propuesto parece contener la fisica esencial y necesaria para 
resultar iitil como un modelo prototipico destinado a describir fenomenos 
anisotropicos observados en la naturaleza. 

• Este modelo es mas conveniente para propositos simulacionales en flui- 
dos con arrastre que otros, como, por ejemplo, simulaciones estandar de 
dinamica molecular, que son mas costosas. 

• El modelo es la extension natural a casos anisotropicos lejos del equilib- 
rio del conocido fluido de Lennard-Jones, que describe acertadamente el 
comportamiento de gases reales como el Ar. En efecto, nuestro modelo 
se reduce a este para campo cero. De otro modo, exhibe una corriente 
neta de particulas y estructuras macroscopicas franjeadas por debajo de 
un punto critico. 

• En cuanto a las propiedades cineticas, y a pesar de la anisotropia inher- 
ente del modelo en las etapas tardias de la descomposicion espinodal, la 
nucleacion se produce por coagulacion tipo Smoluchowski y agregacion de 
Ostwald, las cuales son conocidas por describir los procesos del nucleacion 
de fluidos reales en equilibrio. 

• Inesperadamente, encontramos que el comportamiento critico del modelo 
es consistente con la clase de universalidad de Ising (equilibrio), pero no 
con la del DLG. La principal razon de este aparente desacuerdo y de otros 
detallados en este capitulo, podria ser la falta de simetria particula-hueco 
en el modelo Lennard-Jones con arrastre. 
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• De ahi inferimos que la discretizacion espacial puede cambiar de forma 
significativa no solo las propiedades morfologicas y cineticas sino tambien 
las propiedades criticas. Esto contrasta fuertemente con el concepto de 
universalidad en sistemas en equilibrio, donde el comportamiento de estos 
cerca de los puntos cn'ticos es independiente de los detalles microscopicos. 

• Finalmente, creemos que este modelo motivara nuevos experimentos de de- 
scomposicion espinodal anisotropica, ademas de nuevos trabajos teoricos. 
Dc hecho, nuestras observaciones, por ejemplo, en las propiedades es- 
tructurales y en las primeras etapas de la cinetica de segregacion de 
fases, son facilmente accesibles por experimentos de espectroscopfa y mi- 
crocalorimetricos, respectivamentc. 

Muchas cuestiones quedan aiin por desarroUar. Investigaciones adicionales, las 
cuales estan siendo elaboradas actualmente, involucran estudios detallados de 
las propiedades interfaciales y de los efectos de tamaiio finito. Otros problemas 
que seria conveniente abordar en un futuro son dilucidar de forma precisa el pa- 
pel de la simetn'a hueco-partfcula y de la forma del potencial de interaccion en 
las propiedades cn'ticas; analizar las etapas tardias del proceso de segregacion 
de fases, ya estudiadas en otros sistemas en y fuera del equilibrio, |107l 1108) 
y [1031 1109) . respectivamentc; y estudiar, desde el punto de vista de la teorfa 
de vidrios las caracteristicas morfologicas y estructurales de las fases solida y 
vidriatica (solido policristalino) que presenta nuestro modelo. 

Capitulo O En este capftulo proporcionamos la base adecuada para nuestro 
estudio de gases granulares desarrollada en los Capitulos|6]y[71 Discutimos algu- 
nas controversias, problemas abiertos y avances recientes en fluidos granulares, 
en particular, de la aplicabilidad de teorfas hidrodinamicas. Tambien presenta- 
mos los ingredientes basicos de los metodos computaciones implementados en 
los siguientes capftulos. 

Capitulo O Se presentaron calculos analiticos, ademas de simulaciones de 
dinamica molecular por eventos, para un granulado monodisperso confinado en 
un recipiente anular sin gravedad y forzado por una pared "termalizada" . Con- 
sideramos el limite cuasielastico en una descripcion hidrodinamica que emplea 
coeficientes de transporte tipo Chapman-Enskog [I35j, derivables de primeros 
principios. DesarroUamos un estudio exhaustivo de las inestabilidades de agre- 
gacion heterogenea de partfculas, ruptura de simetrfa y separacion de fases que 
aparecen en el modelo, y que nos permiten determinar con exactitud su dia- 
grama de fases. Para comprobar y complementar las predicciones teoricas im- 
plementamos simulaciones de dinamica molecular por eventos para este sistema. 
Nuestro estudio nos Ueva a desarrollar algunas consideraciones sobre la ffsica de 
las inestabilidades que tienen lugar el modelo. De este capitulo concluimos: 

• El enfoque hidrostatico quedo completamente determinado por tres para- 
metros adimensionales: el de disipacion de energfa, la fraccion de area 
ocupada por los granos y la razon de aspecto del sistema. 
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• Nuestras simulacioncs dc dinamica molecular mostraron (mas claramcnte 
para un niimero suficientemente elevado de parti'culas) la formacion de 
estados agregados con simetria azimutal alejados de la pared termica 
que fuerza al granulado. A pesar del (rclativamente) escueto mimcro de 
particulas empleado, estos estados anulares fueron descritos de forma pre- 
cisa y para densidades moderadas por la solucion numerica de las ecua- 
ciones hidrostaticas propuestas, las cuales emplean relaciones constitutivas 
tipo Enskog. 

• Mediante un analisis de estabilidad marginal se obtuvieron soluciones 
que bifurcan del estado mas basico (el estado anular). Dicho de otro 
modo, estas soluciones predicen la ruptura espontanea de simetria de 
aquellos estados con simetria azimutal (inestables respecto de pequefias 
perturbaciones), hecho que ocurre solamente en un determinado rango 
de parametros. Determinamos las curvas de estabilidad marginal y las 
comparamos con las fronteras de la region espinodal (de compresibilidad 
negativa) del sistema. 

• Se concluye que el mecanismo fisico de la inestabilidad de separacion de 
fases es la compresibilidad negativa del granulado en la direccion azimutal, 
causada por las perdidas de eiiergi'a en las colisiones. Matematicamente, la 
separacion de fases se manifiesta en la existencia de soluciones adicionales 
en cierta region del espacio de parametros a la ecuacion de densidad. 

• Tambien las simulaciones mostraron la existencia de separacion de fases, 
predichas por la descripcion hidrodinamica, pero enmascaradas por grandes 
fluctuacioncs espacio-temporales. Midiendo la distribucion de probabili- 
dad para la amplitud del modo de fundamental del desarrollo de Fourier 
del espectro azimutal de la densidad de particulas, identificamos con clar- 
idad el cambio de estados simetricos a estados con separacion de fases. 

• Ademas encontramos que la region del espacio de parametros asociada a la 
inestabilidad, predicha por la descripcion hidrodinamica, se localiza dentro 
de la region de separacion de fases observada en nuestras simulaciones. 
Esto implica la presencia de una region binodal (de coexistencia de fases) 
donde el estado anular es metaestable. 

• Esperamos que las inestabilidades dc agregacion, ruptura de simetria y 
separacion de fases scan facilinentc observables en experimentos. De he- 
cho, el modelo aqui propuesto es expcrimentalmente accesible, de modo 
que podrfa motivar inievos experimentos en estas inestabilidades si atra- 
jese la atencion de la comunidad pertinente. 

Es claro que la comparacion dc las prcdicciones de nuestra teoria y simulaciones 
con experimentos seria altamente interesante. Al proponer la geometria anular 
esperamos motivar estudios experimentales de estas inestabilidades que podrian 
resultar ventajosos en esta geometria. Un dispositivo experimental anular no 
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presenta paredes laterales (como en una geometria rectangular) y por tanto evi- 
taria las innecesarias, e inevitables por otro lado, perdidas por colision con las 
partfculas. Un posible experimento podria emplear esferas metalicas deslizando 
sobre una superficie cuasi-plana (ligeramente concava) . La pared termica podrfa 
ser implementada por una anillo ligeramente excentico (y quizas rugoso) rotando 
rapidamente. Igualmente interesante para estudios subsiguientes serfa determi- 
nar la naturaleza de las soluciones que se bifurcan (^sub- o supercriticamente?) 
del estado agregado con simetria azimutal. Otros estudios podn'an explotar las 
similitudes entre la fenomenologi'a observada en nuestro modelo y la de los anil- 
los planetarios, en los cuales, tambien se observan inestabilidades de agregacion 
y de ruptura de simetria, entre otras. 

Capitulo^ En este capi'tulo tratamos con hidrodinamica granular y fluctua- 
ciones en un sencillo gas granular bidimensional donde, debido a densidades ele- 
vadas (y no a inelasticidad elevada) , descripciones continuas basadas en primeros 
principios dejan de tener validez. Los enfoques hidrodinamicos fueron puestos 
a prueba en un caso extremo empleando simulaciones de dinamica molecular 
por eventos: cuando se forman agregados granulares cuya densidad es proxima 
al valor de empaquetamiento maximo (denominados macropartfculas). El sis- 
tema modelo considerado esta formado por discos rfgidos inelasticos encerrados 
en un recipiente circular cuya pared esta termalizada forzando al granulado y 
Uevandolo lejos del equilibrio. Las principales conclusiones se resumen en: 

• Simulaciones de dinamica molecular mostraron, a causa del enfriamiento 
por colisiones, un agregado fuertemente empaquetado de partfculas, con 
una interfase cuasicircular abrupta, con fluctuaciones espaciales y ais- 
lado de la pared termica por un gas muy diluido. La formacion de la 
macroparticula requiere un cierto niimero de partfculas y la transicion, en 
funcion de la densidad, desde un estado (cuasi) homogeneo a la macroparti- 
cula es abrupta pero continua. 

• Se resolvio numericamente un conjunto de ecuaciones granulares hidrostati- 
cas que emplean como relaciones constitutivas las fenomenologicas de 
Grossman e< [156) . Encontramos que, en un amplio rango de parametros, 
el ajuste — incluyendo la parte de maximo empaquetamiento — con la 
solucion con simetria azimutal de las ecuaciones hidrostaticas es casi per- 
fecto. Tambien demostramos que, para el mismo sistema, el empleo de 
relaciones constitutivas tipo Enskog [135 , derivables de primeros princip- 
ios, no describe adecuadamente — mediante comparacion directa con los 
resultados de simulacion — la region de mayor densidad del granulado. 

• Estos resultados, que estan de acuerdo con resultados previos en ge- 
ometrfas planas con y sin gravedad, extienden la aplicabilidad de las 
relaciones constitutivas de Grossman et al. a coeficientes de restitucion 
ligeramente menores y a geometrfas no planas. 

• Un analisis de estabilidad marginal demostro — independientemente de las 
relaciones constitutivas empleadas — que no hay soluciones estacionarias 
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que bifurquen linealmente de la de simetria azimutal, esto es, no hay 
estados que manifiesten ruptura de simetn'a alguna. Esto esta de acuerdo 
con la persistencia de macroparticulas circulares y cuasicirculares que se 
observan en nuestras simulaciones. 

• Estudiamos, por primera vez, las fluctuaciones de las macroparticulas, 
midiendo la distribucion de probabilidad radial del centro de masas del 
sistema. Las fluctuaciones resultan ser gaussianas. 

• Las fluctuaciones gaussianas observadas del centro de masas sugieren una 
descripcion mesoscopica tipo Langevin en terminos de una macropartfcula 
confinada un potencial atractivo de naturaleza hidrodinamica y forzada 
por ruido de niimero finito — queremos decir que el agregado efectiia un 
movimiento browniano dentro del sistema. Se determine la correspondi- 
ente ecuacion de Langevin. 

• En contra de lo que ocurre en sistemas en equilibrio — y en muchos lejos 
del equilibrio — , en los cuales la magnitud relativa de las fluctuaciones 
decrece al incrementar el mimero de partfculaj^. encontramos que las 
fluctuaciones persisten al incrementar el niimero de partfculas. 

Aiin quedan algunos problemas en nuestro modelo circular. Una cuestion que 
deberfa ser tratada es si existen inestabilidades adicionales en este sistema mod- 
elo. Algo de intuicion podria obtenerse, por ejemplo, mediante medidas de or- 
dcn posicional y orientacional del agregado, siguiendo la analogfa con la teorfa 
KTHNY '119] para fusion o licuefaccion en sistemas en equilibrio. Concreta- 
mente, estas podrfan Uevarse a cabo monitorizando las respectivas funciones de 
correlacion. Es posible que tanto una solucion general dependiente del tiempo 
como un estudio de estabilidad no-lineal contribuyan a la biisqueda. Otros 
posibles trabajos futuros son el estudio de inestabilidades en modelos tridimen- 
sionales relacionados (en el Apendice |D] se incluye una primera aproximacion a 
este tema) y de la fenomenologfa mostrada por gases granulares polidispersos 
en este tipo de sistemas. 

A pesar de que comiinmente los apendices se relegan a un segundo piano, 
creemos que nuestros Apendices [Bl v [Dl merecen sus propias conclusiones. 

Apendice\B[ Recordamos que en este apendice presentamos un modelo difu- 
sivo y con arrastre cuyo proposito es comprender mejor el balance entre campo 
externo y baiio termico que ocurre en cl DLG y en otros modelos relaciona- 
dos (concretamente, los modelos reticulares estudiados en el Capftulo [3]). Los 
principales resultados son: 

• Estudiamos la estabilidad de la interfase (orientada segiin la direccion del 
campo) en el li'mite de campos saturantes, encontrando un cambio de fase 



^*Para un gas ideal en equilibrio: ff ~ OfN '^^], donde a es la magnitud relativa de las 
fluctuaciones y N es el numero de parti'culas. 
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de segundo orden geometrico o topologico en funcion de la conectividad de 
los sitios del reticulo. Para bajas conectividades existe orden, dominan los 
efectos termicos y la interfase es estable, mientras que al incrementar la 
conectividad (por encima de un valor critico) el orden desaparece, domina 
el campo externo y la interfase se vuelve inestable. 

• Este cambio de fase geometrico confirma que el balance campo-temperatura 
no es tan sencillo como la dinamica microscopica sugiere, es decir, no es 
una mera superposicion lineal de ambos. Este hecho es una manifestacion 
clara de una dinamica colectiva emergente. 

• Ademas se determinaron los parametros criticos del cambio de fase. 

Apendice O En este apendice desarroUamos una serie de estudios teoricos 
para mostrar la existencia de inestabilidades de agregacion y ruptura de simetria 
en flujos granulares tridimensionales. De hecho, se trata de una sencilla ex- 
tension a tres dimensiones de la geometn'a circular estudiada en el Capi'tulo [7l 
Asf, el modelo es el de un gas de esferas n'gidas colisionando inelasticamente en 
el interior de un cilindro cuya pared curva esta termalizada. Para este estudio, 
se hizo necesario proponer un nuevo conjunto de relaciones constitutivas. 

• Derivamos relaciones constitutivas para gases granulares tridimensionales 
empleando argumentos de volumen lihre en la vecindad de la densidad de 
maximo empaquetamiento, y sugerimos una interpolacion entre el Ifmite 
de maxima densidad y las relaciones usuales en el li'mite diluido. Para ello 
nos basamos en la estrategia seguida por Grossman et al. [156j en dos 
dimensiones, que emplea tecnicas de volumen excluido para determinar 
las relaciones constitutivas. 

• Mediante los metodos desarroUados en los Capi'tuloslHly Elpronosticamos 
y caracterizamos las inestabilidades de agregacion y ruptura de simetria 
en un recipiente cih'ndrico. La presencia de esta ultima supone una gran 
diferencia con lo observado en el Capi'tulo [71 donde su homologo bidimen- 
sional no presentaba ruptura de simetria alguna. 
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